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ABSTRACT 

This  thesis  uses  the  state  transition  technique  to 
obtain  the  z-transforms  of  digital  controllers  which,  when 
inserted  into  the  particular  system  under  study,  provide 
an  output  response  that  is  deadbeat.  Various  methods  of 
graphically  representing  the  system  including  the  digital 
controller,  writing  the  state  transition  equations  from  the 
flow  graphs  and  solving  the  equations  are  described. 

A  digital  computer  program  was  written  to  solve  the 
state  equations  for  third- and  fourth-order  systems.  The 
program  was  used  to  obtain  the  transforms  of  digital  control¬ 
lers  for  a  number  of  systems  and  the  results  confirmed  on  the 
PACE  231-R  analog  computer. 

A  method,  proposed  by  Y.  J.  Kingma  and  programmed  in 
FORTRAN  II  by  D.  H.  Kelly,  of  solving  for  the  roots  of  a 
polynomial  using  Routin’ s  stability  criterion  was  rewritten 
into  FORTRAN  IV  and  double  precisioned  to  provide  results 
with  greater  accuracy.  The  program  will  solve  ill-conditioned 
polynomials  and  will  provide  results  as  accurate  as  those 
given  by  the  SHARE  Library  program  C2-3332  based  on  Newton 
and  Muller  approximations. 
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CHAPTER  I 


INTRODUCTION 

In  the  synthesis  and  design  of  feedback  control  systems 
the  output  response  must  satisfy  some  predetermined  performance 
Some  of  the  criteria  used  in  the  design  of  both  sampled-data 
and  all-continuous  systems  are  maximum  overshoot,  the  time 
for  the  error  to  reach  its  first  zero,  the  time  to  reach 
maximum  overshoot,  settling  time  and  the  frequency  of  oscil¬ 
lation  of  the  transient.  If  the  design  requirements  are 
such  that  a  rigid  performance  specification  must  be  met, 
the  usual  methods  of  compensation  become  cumbersome,  time- 
consuming  and  depend  on  the  experience  of  the  designer.  This 
is  especially  true  when  working  with  sampled-data  systems. 

The  compensation  of  sampled-data  systems  to  obtain 
the  desired  performance  can  be  accomplished  by  insertion  of 
continuous-data  networks  in  the  forward  or  feedback  paths 
of  the  system.  Since  the  network  is  connected  to  the  con¬ 
trolled  process,  the  z-transform  of  the  product  of  the  two 
transfer  functions  is  a  single  function.  Analytically,  this 
makes  it  difficult  to  determine  the  effect  of  various  networks 
on  the  overall  performance.  As  an  alternative,  digital 
compensation  can  be  implemented  through  the  use  of  digital 
controllers.  Digital  controllers  are  in  general  more 
versatile  than  continuous-data  networks  and  better  system 
performance  can  thus  be  achieved.  They  can  be  used  in 
continuous-data  systems  as  well  as  sampled-data  systems. 

Digital  controllers  may  be  synthesized  using  passive 
networks  preceded  and  followed  by  sampling  switches,  by 
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digital  computer  techniques  or  by  employing  both  analog  and 
digital  components.  Its  main  function  is  to  process  a  sequence 
of  numbers  equally  spaced  in  time  into  a  command  signal 
which  is  then  applied  to  the  control  process. 

Using  this  method  of  compensation,  quite  stringent 
performance  specifications  can  be  met,  such  as  minimal  and 
deadbeat  response.  Minimal  response  requires  the  system  to 
have : 

1.  Zero  steady-state  error  at  the  sampling  instants 
for  a  specified  input  signal. 

2.  A  rise  time  equal  to  a  minimum  number  of  sampling 
periods . 

3.  A  finite  settling  time  measured  at  the  sampling 
instants . 

The  design  requirements  for  deadbeat  response  are  more  rigid 
than  for  minimal  response.  They  are  as  follows: 

1.  The  system  must  have  zero  steady-state  error  after 
a  minimum  number  of  sampling  instants  with  no 
intersampling  ripple  for  a  specified  input  signal. 

2.  The  transient  response  should  be  as  fast  as  possible, 
with  a  finite  settling  time  measured  at  the  sampling 
instants . 

Deadbeat  response  is  a  form  of  minimal  response  with  the 
further  restriction  of  no  intersampling  ripple.  The  ensuing 
chapters  will  demonstrate  the  synthesis  procedures  for 
obtaining  deadbeat  response  using  digital  controllers. 

To  design  a  digital  controller  for  deadbeat  response, 
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the  state  transition  approach  will  be  used  as  it  provides  a 
simple  and  systematic  method  for  continuous-data  as  well  as 
sampled-data  systems.  The  procedure  will  be  developed  using 
the  transition  flow  graphs  and  equations,  and  a  digital 
program  will  be  described  for  solving  the  state  equations. 

The  method  is  easy  to  use  and  requires  little  design 
experience.  Although  the  method  can  be  applied  to  systems 
having  real,  finite  poles  with  assurance  of  obtaining  the 
desired  results,  a  judicious  choice  of  sampling  rate  is 
required  to  obtain  a  stable  output  when  the  system  contains 
complex  conjugate  roots. 

Since  a  knowledge  of  the  roots  of  the  plant  polynomial 
is  sometimes  required,  a  method  of  determining  them  using 
Routh's  stability  criterion  is  described.  Knowing  the  roots, 
the  designer  can  then  choose  an  appropriate  sampling  rate 
to  obtain  the  desired  deadbeat  response. 


i  £.  w  h  8  b;  e  -  8  jo  Li  ill  ^  no  o  iol  boritsm  ol  Jamecte^e  fens  slqmte 

f  09  -.  I  •'  *vIos  rr  .  I1  f  '  1 1  .  r;  OK 

ns  Is  9  b  ■''  C  3:;  :Jtc/po  bn 3  s  j  .  ■  ol  >o  d->m  eriT 

m  •.■'■■ 

srid  snlnisddo  lo  son  tes  ridiw  se C oq  ^Inll  *Is9tt 

' 

-  8  £  '  L 

enlsdnoo  made^e  ;  j  nariw  dt/q;tuo  [<  >o  •  s  nl  ido  od  upei 

§ni8u  msrid  gnlnlfitied-  b  lo  f  v  d  a.  £*'?$•  *i  tupsn:  armlet  .nos  el 


. 


CHAPTER  II 


Z-TRANSFORMS  OF  DIGITAL  CONTROLLERS  USING  THE  STATE  TRANSITION 
TECHNIQUE 

There  are  two  major  approaches  to  system  analysis  and 
synthesis:  l)  the  use  of  block  diagrams  and  transfer  functions, 

and  2)  the  state  transition  or  state  variable  approach  using 
signal  flow  graphs.  The  state  transition  approach  will  be 
used  in  this  thesis  since  the  system  can  be  represented  by  a 
set  of  first-order  differential  or  difference  equations 
describing  the  state  variables  of  the  system.  The  continuous 
elements  are  described  by  differential  equations,  the  discrete 
elements  by  difference  equations,  and  the  sample-and-hold 
elements  by  discrete  state  transition  equations. 

Once  mathematical  expressions  have  been  obtained  for 
all  the  system  elements,  the  state  transition  flow  graph  can 
be  drawn.  With  the  aid  of  the  flow  graph,  the  state  transition 
equations  can  readily  be  written.  Solving  these  equations 
manually  is  time  consuming,  but  since  they  are  linear  differ¬ 
ence  equations,  they  can  be  solved  using  a  digital  computer. 

2.1  BLOCK  DIAGRAM  REPRESENTATION  OF  CONTROL  SYSTEM 

If  a  digital  controller  is  used  to  compensate  a 
sampled-data  system,  its  block  diagram  representation  consists 
of  the  Laplace  transform  of  the  controller  followed  by  a 
sampling  switch  synchronized  with  the  other  system  switches. 

The  block  diagram  is  drawn  as  shown  in  Figure  2-1. 

If  the  system  to  be  compensated  is  continuous,  the 
block  diagram  representation  of  the  digital  controller 
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FIGURE  2-1.  SAMPLED-DATA  SYSTEM 


Digital  Controller 
I - I 


FIGURE  2-2.  CONTINUOUS-DATA  SYSTEM 
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consists  of  the  Laplace  transform  of  the  controller  preceded 
and  followed  by  synchronous  switches  and  followed  by  a  zero- 
order  hold  as  shown  in  Figure  2-2. 

2.2  STATE  VARIABLES  AND  STATE  TRANSITION  EQUATIONS  OF 
CONTINUOUS  ELEMENTS 

An  n-th  order,  linear,  continuous  system  can  be  repre¬ 
sented  by  a  set  of  first-order  differential  equations  of  the 

form:  dx^  =  f^(xj,  X2*...xn,  r^,  r2....rm)  where  the  x's  are 

dt 

the  set  of  n  system  variables  (state  variables);  the  r's  are 
the  m  system  inputs;  the  f's  are  the  functional  relations. 

In  vector-matrix  notation  these  equations  can  be 
written  as 

dx  =  Ax  f  Br  . (2-1 ) 

dt 

where  x  is  a  column  matrix  of  the  state  variables  and  r  is  a 
column  matrix  of  the  inputs.  The  equation  describes  the 
transition  of  state  of  the  system's  state  variables  from 
x^(t)  to  Xj_(t  )f  dXj_(t )  when  the  time  is  changed  from  t  to 
t-t-dt.  The  transition  is  continuous  since  dt  is  an  infini- 
tesmal  time  change. 

Obtaining  the  Laplace  transform  of  (2-1)  and 
rearranging  the  result  yields 

X(s)  -  (si- A)  ^(tj)  f  ( sI-A )  ^BR ( s ) 
where  x(tQ)  is  a  column  matrix  of  the  initial  conditions  of 
x(t)  evaluated  at  t  =  tj.  The  inverse  Laplace  transform  of 
X(s)  yields  ^ 

x(t)  -  J(t-t0)x(tQ)  +  |  J(t-r)Br(T)dr 

-40 


(2-2) 
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where  J(t-t0)  =  =  J^^((sT-A)  ^")  and  the  second 

term  on  the  right-hand  side  of  (2-2)  is  the  convolution 
integral.  This  set  of  equations  is  defined  as  the  state 
transition  equations  and  J(t-tQ)  is  referred  to  as  the  state 
transition  matrix  of  the  dynamic  system. 

Using  the  vector-matrix  method  provides  simplified 
notations  for  obtaining  a  solution,,  but  the  signal  flow  graph 
approach  provides  a  more  straightforward  method  of  writing 
the  state  transition  equations.  The  state  transition  flow 
graph  of  a  system  is  the  analog  computer  simulation  flow  graph 
with  the  initial  conditions  applied  as  input  signals  at  the 
corresponding  output  nodes  of  the  integrators.  Using  the 
flow  graph  and  Mason's  gain  formula^  the  state  transition 
equations  can  readily  be  written. 

The  analog  simulation  flow  graph  of  the  continuous 
elements  may  be  programmed  by  three  different  methods.  The 
two  most  common  methods  are  direct  programming  and  iterative 
programming,  with  the  third  method  being  parallel  programming. 
The  method  one  should  use  is  decided  by  the  form  of  the  given 
transfer  function.  If  the  transfer  function  is  not  in  fact¬ 
ored  form,  direct  programming  is  used,  but  if  it  is  in  fact¬ 
ored  form  or  can  be  factored,  then  iterative  programming  is 
the  easiest  method  to  use.  Parallel  programming  is  used  if 
the  transfer  function  can  be  expressed  in  partial  fraction 
form.  Using  iterative  programming  reduces  the  number  of 
terms  in  the  state  transition  equations  and  thus  the  number 
of  calculations  required  is  reduced.  If  the  calculations 
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are  to  be  carried  out  manually,  this  method  not  only  saves 
time,  but  reduces  the  chances  of  incurring  errors. 

To  illustrate  direct  and  iterative  programming,  an 
example  will  be  developed  using  both  methods.  The  chosen 
plant  has  the  transfer  function  G(s)  =  l/( s3  +  3s2+2s )  or  in 
factored  form  G(s)  =  l/s  ( s  +  1 ) ( Sf2) . 

Using  the  direct  method,  the  expression  for  G(s)  can 
be  written  as 

(s3  +  3s2+2s  )C(  s)  =  R(s)  . (2-3) 

where  C(s)  and  R(s)  represent  the  output  and  input  of  the 
plant  respectively.  Choosing  the  state  variables  Xj_(s)  as 

X1(s)  =  C ( s ) 

sX^s)  =  x2(s) 

S2Xl(s)  =  X3(s) 

S3x1(s)  =  X2t(s) 

and  substituting  into  expression  (2-3)  yields 

X4(s)  =  R(s)  -  3X3(s)  -  2X2(s). 

To  draw  the  analog  simulation  flow  graph,  the 
integrators  are  represented  by  s--*-  and  X j,  Xg,  X3  are  chosen 
as  the  nodes.  The  flow  graph  is  shown  in  Figure  (2-3). 


FIGURE  2-3.  FLOW  GRAPH  USING  DIRECT  METHOD 
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Using  the  iterative  method,  the  transfer  function  is 


separated  into  its  factors  and 
for  each  factor  thusly 

Sf  2 
1 


s  +  1 


a  separate  flow  graph  is  drawn 


To  draw  the  complete  flow  graph,  the  separate  flow 
graphs  are  joined  together  using  unity  gain  branches  as  shown 
in  Figure  (2-4)  with  the  output  nodes  of  each  graph  labelled 
as  Xj,  and  X^. 


FIGURE  2-4.  FLOW  GRAPH  USING  ITERATIVE  METHOD 

The  state  transition  flow  graph  can  be  drawn  directly 
from  the  analog  simulation  flow  graph.  The  output  nodes 
of  the  integrators  are  chosen  as  the  state  variables  and  the 
initial  conditions  are  applied  as  inputs  to  the  output  nodes. 

2.3  SAMPLE-AND-HOLD  ELEMENTS 

The  sample-and-hold  elements,  the  sampling  switch  and 
the  zero-order  hold,  are  normally  positioned  between  the 
digital  controller  and  the  continuous  elements.  The  output 
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of  the  digital  controller  is  sampled  by  a  switch  closed  for  a 
short  duration,  short  enough  that  the  input  to  the  zero-order 
hold  appears  as  an  impulse.  The  zero-order  hold  retains  the 
value  until  the  switch  closes  at  the  next  sampling  instant 
T  seconds  later.  Thus  the  output  of  the  zero-order  hold,  the 
input  to  the  continuous  elements,  is  a  step  function. 

If  the  input  and  output  are  denoted  by  m(t)  and  h(t) 
respectively,  then  h(kT+)  *  m(kT)  where  kTf  denotes  the  k-th 
sampling  period  from  time  kT  to  kTt-'C  for  0<rC'^T.  The  Laplace 
transform  of  a  step  function  is  l/s,  thus  on  the  state 
transition  flow  graph  the  sample-and-hold  elements  are  repre¬ 
sented  by  s"1. 

2.4  VARIABLE  GAIN  CONCEPT  OF  DIGITAL  CONTROLLERS 

To  complete  the  drawing  of  the  state  transition  flow 
graph,  the  digital  controller  must  be  described  by  a  mathema¬ 
tical  expression.  The  variable  gain  concept  will  be  employed 
for  this  purpose. 

Referring  to  Figure  2-1,  the  output  of  the  zero-order 
hold  is  h(t)  and  the  controller  actuating  signal  is  e(t). 

The  digital  controller  output  is  a  train  of  impulses  with 
values  at  the  sampling  instants  which  are  equal  to  the  zero- 
order  hold  output  at  the  sampling  instants.  Under  this 
condition  the  z-transfer  function  can  be  written  as 

D(z)  =  H(z)/E(z) 

where  H(z)  =  h(0+)  +  h(Tf)z"1t  h(2Tf)z_2f - f  h(nT+)z"n  and 

E(z)  =■  e(0+)  t  e(T*)z-1  +  e  (2T+ )  z_2+- .  .  .  .  +  e(nT+)z~n. 
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Using  the  state  variable  approach,  the  digital  controller 
may  be  replaced  by  a  variable  gain  amplifier  with  gain  K(kT). 
The  amplifier  will  have  different  values  during  different 
sampling  periods  since  during  any  one  sampling  period  only 
one  term  in  each  of  the  equations  for  H(z)  and  E(z)  is  of 
interest.  The  gain  during  the  first  sampling  period  is  given 
by  K(0)  =  h(0+ )/e (0^ ) ,  the  second  sampling  period  by 
K(T)  =  h( 1T+ )/e ( 1T+ )  and  so  on.  Therefore  the  gain  during 
the  (k+l)-th  sampling  period  is  given  by  K(kT)  =  h(kT+' )/e (kTf )  . 
Thus,  when  drawing  the  state  transition  flow  graph  for  the 
time  interval  kT  <t  ^.(k+-l)T,  the  digital  controller  will  be 
represented  by  the  variable  gain  amplifier  with  gain  K(kT). 

2.5  STATE  TRANSITION  FLOW  GRAPH 

With  the  preceding  element  descriptions,  the  complete 
state  transition  flow  graph  may  be  drawn.  As  previously 
mentioned,  to  represent  the  continuous  elements,  the  analog 
simulation  flow  graph  is  modified  by  adding  initial  conditions 
to  the  output  nodes  of  the  integrators.  The  initial  condi¬ 
tions  are  represented  by  a  unity  gain  branch  from  x(tQ)  to 
x(tj)  and  from  x(tj)  to  the  state  variable  X(s)  by  a  branch 
with  gain  l/s.  The  reason  for  this  is  that  the  initial 
conditions  appear  initially  as  step  functions  at  the  outputs 
of  the  integrators. 

To  illustrate  the  final  flow  graph  an  example  will  be 
described  using  the  closed-loop  system  of  Figure  2-5.  The 
analog  simulation  flow  graph  of  G(s)  is  shown  in  Figure  2-4. 
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FIGURE  2-5.  CLOSED -LOOP  SYSTEM 

Combining  this  with  the  variable  gain  concept  of  the  digital 
controller  and  the  sample-and-hold  element  representation, 
the  complete  state  transition  flow  graph  may  be  drawn  as 
shown  in  Figure  2-6.  The  digital  controller  is  described 
by  K(t0) . 


2.6  STATE  TRANSITION  EQUATIONS 

The  state  transition  equations  can  be  written  directly 

•3 

from  the  flow  graph  using  Mason's  gain  formula  .  The 
equations  are: 


1 
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x1(s) 


’  x1(to)[1/s-K(t0)s  k/& 


1 
4  X 


4  X, 


(t0)fs'2(l43s"1+2s'2)/A 


(t0)[s  3 /A  4  r(t0)K(t0)|js  ^  /A 


X2(s)  -  x1(t0)[^-K(t0)s  3/^J  +  x2(tQ)  s  1(l-2s  1)/& 


1 
4  X 


(t0)[s-2/A] 

[- 


s'3/A 


4  r(t0)K(t0) 

X,(s)  =  x1(tQ)| -K(t0)s'2/A  4  x3(t0)|s-1/A 


. (2-4) 


4  r(t0)K(t0) 


s'2/A 


where  A-  (l+3s~"*'42s’2) . 

The  digital  controller  starts  operating  at  time  t0-0 
seconds  and  the  output  of  the  plant  is  continuous  until  the 
next  sampling  period,  T  seconds  later.  The  values  of  the 
state  variables  at  time  T  become  the  initial  conditions  for 
the  next  sampling  period.  Taking  the  inverse  Laplace  trans¬ 
form  of  equations  (2-4)  and  replacing  tQ  by  kT,  t  by  (k4l)T 
and  K(kT)  by  Kk  the  equations  become 


x 


1 


(k4l)T 


-Kk(T/2-3/44e  T-e“2T/4)  +  1  x2(kT) 


4  0 . 5 ( l-e  2T^x2(kT)  ¥  0. 5( 1+e  "x-2e  x)x3(kT) 


-2T  ^  -T 


+  r(kT)Kk(T/2-3/4+e"T-e“2T/^) 


(2-5) 


x, 


(k4l)T 


=  -0.5Kk(lte"2T-2e"T)x1(kT)  +  (e  2T)x2(kT) 


4  (e"T-e"2T)x3(kT)  4  0. 5Kk( lte“2T-2e'T)r (kT) 


-2T  o.-T 


x- 


(k+l)T 


-  -Kk(l-e"T)x1(kT)  4  (e_T)x3(kT)  4  ( l-e"'1  )Kkr (kT) 


-T 


-T’ 


For  deadbeat  response,  the  system  error  must  be  zero 
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for  t  >  nT  where  n  is  the  smallest  possible  positive  integer. 
For  a  third-order  system,  three  is  the  smallest  value  n  can 
have.  The  reason  for  this  is  that  there  are  three  conditions 
which  must  be  satisfied,  the  values  of  the  three  state  vari¬ 
ables.  The  state  variables  are  functions  of  the  variable 
gain  where  k  takes  on  the  values  0,  1  and  2.  Letting  k 
take  on  these  values  produces  three  equations  in  three 
unknowns.  Applying  a  unit  step  input  to  the  above  system,  the 
output  of  the  last  integrator  (x-^)  would  be  unity  for  the 
system  to  have  zero  steady  state  error.  In  order  that  this 
output  does  not  change,  the  remaining  state  variables  (x2  and 
x^)  must  be  zero  at  the  end  of  the  third  sampling  period. 

Proceeding  with  the  solution  for  a  unit  step  input 
and  a  sampling  period  of  T=  1  second,  the  equations  become 


x-jftkt-ljTj-  (-.0840461^  +  l)x1(kT)  +  (,432332)x2(kT) 

+  (.  199786  )x3(kT)  +■  (.084046)Kk 


(k>l)T 


(-.199786Kk)x1(kT)  4-  (.135335)x2(kT) 
4.  (,23254)x3(kT)  +  (.199786)1^ 


x3j(k+l)Tl-  (-.63212Kk)x1(kT)  +  (.  36788  )x3(kT) 

+  (,63212)Kk 


(2-6) 


Letting  k=0  in  (2-6)  yields  the  values  of  the  state 
variables  at  the  end  of  the  first  sampling  period  which  are 

x,(T)  =  .084046KQ 
x2(T)  =  ,199786Kb 
x3(T)  =  .63212Kq 


(2-7) 
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Setting  k=l  in  (2-6)  and  substituting  (2-7)  into  (2-6) 
results  in  the  equations 

x1(2T)  =  -.00706373^1^  +  .084046K1  +  .296709K0 

x2(2T)  =  -.0167912K0K1  +  .  199786K-L  +  .174031K0  . (2-8) 

x3(2T)  «  -.05312721^1^  +  .  63212K2  +  .  2325^  Kq 

Making  a  final  substitution  of  k=2  in  the  equations 
(2-6)  and  substituting  the  expressions  (2-8)  into  (2-6) 
produces  the  required  equations  which  are 

xa(3T)  =  .000593678%%%  -  .00706373X2%  -  .0249372%% 

-.0249372%%+  .084046%+  .296709%+  .418407% _ (2-9) 

x2(3T)  =  .00141123%%%  -  .0167912%%  -  .0592783%% 

-  .0146266%%  +  .199786%  +  .174031%  +  .0776283%'' 

x3(3T)  =  .00446513%%%  -  .0531272%%  -  .187556%% 

-.0195444%%  +  .232544%  +  .63212%  +  .0855483%'"' 

For  deadbeat  response  x,(3T)  =  1,  Xg(3T)  =  0,  and 
x3(3T)  =  0.  There  are  now  three  simultaneous  equations  in 
the  three  unknowns  Kq,  and  Kg.  These  equations  are  solved 
by  eliminating  the  unknowns  in  a  particular  sequence,  starting 
with  the  gain  having  the  largest  subscript.  That  is,  multi¬ 
plying  by  the  coefficients  of  the  Kg  terms  in  a  third-order 
system  and  subtracting  the  result,  eliminates  all  terms 
containing  Kg  from  the  equations.  Continuing,  choosing  the 
coefficients  of  the  terms,  multiplying  and  subtracting, 
eliminates  all  terms  with  the  factor  Kj,  leaving  one  equation 
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in  the  remaining  unknown  Kq.  Beginning  the  elimination  of 
unknowns  with  the  coefficients  of  a  term  other  than  %  will 
not  eliminate  all  the  product  terms  containing  that  factor. 

For  example,  multiplying  (2-10)  by  .63212  and  (2-11) 
by  0.199786  and  subtracting  the  result  yields 

-.00534107%%  +  .063549%  +  .031979%  =  0  . (2-12) 

in  which  all  terms  containing  the  factor  %  have  been  elimin¬ 
ated.  Continuing  the  solution  in  this  manner  results  in  the 
required  variable  gains 

%  =  3.659199 

%  =  -2.659186  . (2-13) 

k2  -  2.638810 

Having  obtained  the  variable  gains  the  next  step. is 
to  determine  the  transfer  function  of  the  digital  controller. 

2.7  TRANSFER  FUNCTION  OF  DIGITAL  CONTROLLER 

To  determine  the  transfer  function,  the  error  functions, 
that  is,  inputs  to  the  digital  controller,  must  be  calculated 
at  the  sampling  instants.  The  input  r(kT)  is  a  unit  step  and 
the  output  can  be  determined  from  the  state  variable  equations 
(2-7)  and  (2-8)  for  x^(T)  and  x^(2T)  respectively.  The  output 
is  considered  to  be  initially  zero,  therefore  the  error 
functions  are 

e(0+)  =  r(0)  -  x  (0)  =  1 
e (T* )  =  r(T)  -  x1 (t)  =  .692459 
e(2T+)  =  r(2T)  -  x1(2T)  *  .069043 
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The  outputs  of  the  digital  controller  are  found  from 
the  variable  gain  equation  -  h(kT+ )/e (kT+ )  and  are 
calculated  to  be 

h(0+)  -  KQe(0+)  -  3.659199 
h(T+)  =  K^e (T** )  -  -1.841377 
h(2T+)  =  K2e(2T+)  =  0.182191 

The  complete  transfer  function  can  now  be  written 

D(z)  =  3-659199  -  1.84l377z~*  +  0,l82191z~^ 

1  +  0.692il59z'1  +  0.069043z'2' 

The  z-transform  is  used  since  the' complete  system  is 
a  sampled-data  system  and  only  one  term  is  valid  during  each 
sampling  period. 

The  digital  controller  can  be  realized  on  the  analog 
computer  with  the  aid  of  time  delays  or  on  a  digital  computer 
using  data  storage^. 

As  can  be  seen  from  the  above  example,  the  numbers 
involved  in  the  solution  are  cumbersome  when  solving  the 
equations  on  a  calculator  and  the  possibilities  of  incurring 
errors  are  numerous.  To  obtain  an  accurate  answer  quickly 
a  digital  computer  is  required.  The  next  chapter  is  devoted 
to  a  FORTRAN  program  for  use  on  the  IBM  7040  computer. 
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CHAPTER  III 


PROGRAM  FOR  SOLVING  STATE  EQUATIONS 

The  state  transition  approach  to  calculating  the 
transfer  function  of  a  digital  controller,  which  will  produce 
deadbeat  response  when  used  as  a  compensator,  lends  itself  to 
solution  on  a  digital  computer.  All  calculations  are  step- 
by-step  multiplications,  subtractions  and  divisions  as  can  be 
seen  from  the  example  in  Chapter  II.  The  program  to  be 
described  in  this  chapter  is  not  a  general  program,  but  one 
written  to  solve  specific  types  of  systems. 

3.1  SYSTEMS  SOLVABLE  USING  DIGITAL  PROGRAM 

The  program  was  designed  for  use  with  linear  third- 
and  fourth-order  systems.  The  transfer  functions  of  these 
systems  must  not  contain  any  zeros  and  must  have  at  least 
one  pole  at  the  origin,  that  is,  a  free  integrator.  A 
further  restriction  is  that  the  excitation  function  be  a 
step  function.  For  calculation  purposes,  the  value  of  the 
step  input  can  be  assumed  to  be  unity  since  the  system  is 
linear  and  continuous  between  sampling  instants. 

The  system  flow  graph  must  be  drawn  using  the  itera¬ 
tive  or  direct  methods.  If  the  partial-fraction  method  of 
drawing  the  flow  graph  were  used  then  the  program  would  not 
provide  a  solution.  The  reason  will  be  explained  in  a  later 
section . 

3.2  DIGITAL  PROGRAM  FOR  USE  ON  THE  IBM  70^0  COMPUTER 


The  program,  written  in  FORTRAN  IV  computing  language, 
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performs  each  calculation  carried  out  in  the  solution  of  the 
example  described  in  Chapter  II. 

The  coefficients  of  the  state  variables  are  placed 
in  an  E  or  F  matrix  depending  on  whether  the  system  is  of 
third-  or  fourth-order  respectively.  Using  the  example  of 
Chapter  II,  the  coefficients  of  equations  (2-5)  would  be 
placed  in  the  first  three  rows  of  the  E  matrix.  Making  the 
appropriate  substitutions  and  calculations,  the  values  of 
the  state  variables  in  terms  of  the  variable  gains  are  obtained 
at  the  sampling  instants.  The  coefficients  of  the  variable 
gains  and  products  of  variable  gains,  as  shown  in  equations 
(2-7)  to  (2-11),  are  stored  in  succeeding  rows  of  the  matrix. 
The  final  three  rows  of  the  matrix  contain  the  coefficients 
of  equations  (2-9)  to  (2-11)  which  are  the  three  equations 
in  the  three  unknowns  (the  variable  gains)  that  have  to  be 
solved.  To  solve  these  three  equations,  their  coefficients 
are  first  transferred  to  an  A  matrix,  to  a  B  matrix  if  the 
system  is  of  fourth-order,  where  the  unknowns  are  eliminated 
and  the  coefficients  of  the  reduced  equations  are  stored. 

The  final  row  of  the  matrix  A  contains  the  coefficients 
of  the  equation  in  the  unknown  K^.  The  variable  gains  are 
then  solved  for  along  with  the  values  of  the  errors  at  the 
sampling  instants  and  the  values  of  the  inputs  to  the  plant. 

The  program  shown  below  contains  a  number  of  state¬ 
ments  precededed  by  a  "C",  which  indicates  the  statement  is  a 
comment.  Comment  statements  do  not  affect  the  execution  of 
the  program,  and  are  used  to  indicate  what  calculations  are 


to  follow. 
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Preceding  the  program  is  a  list  of  the  terms  used  in 
the  program. 

CKO,  CK1,  CK2,  CK3 - Variable  gains  Kq,  K-,,  Kg, 


CKN2  - Numerator  of  expression  for  CK2 

CKD2  - Denominator  of  expression  for  CK2 

ERO - error  at  time  t=0* 

ER1 - error  at  time  t=T+ 

ER2 - error  at  time  t=2T+ 

HO - input  to  plant  at  time  t=0+ 

Hi - - - input  to  plant  at  time  t=T+ 

H2 - input  to  plant  at  time  t=2T+ 

T - sampling  period  in  seconds 

N - order  of  plant 


The  flow  chart  of  the  digital  program  shown  below 
is  shown  in  Appendix  A. 

$IBFTC  EQNS 

C  SOLUTION  OF  STATE  EQUATIONS  TO  OBTAIN  THE  EXPRESSION 
C  FOR  A  DIGITAL  CONTROLLER 

COMMON  A(12,8),  B(l3,l6),  E(l2,7),  F(17,15) 

READ  (5,1)  N,T 

1  FORMAT  (I5,E12.3) 

WRITE  (6,2)  N,T 

2  FORMAT ( IX, 15, E10 . 3 ) 

CALL  ARRAYS  (N,T) 

IF  (N.EQ.3)  GO  TO  Hg 

C  OBTAIN  COEFFICIENTS  OF  FOURTH- ORDER  STATE  VARIABLE  EQUATIONS 
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F(5,l)  -  F(l,6) 
DO  3  J  -  2,N 
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3  F(5,J)  -  f(  j,5) 
do4  I=1,B 

it  F(l+5,l)  =  F(I,1)*F(5,1) 

F(6,2)  =  F(l,6) 

DO  5  I  -  2, N 

5  F(l+5,2)  =  F(I,5) 

f(6,3)=f(i,2)*f(5,1)+f(1,3)*f(5,2)+f(i,4)*f(5,3)+f(1,5)* 

1F(5,4) 

DO  7  1-2,  N 

7  F(I+5,3)=F(I,2)*F(5,2)+F(I,3)*F(5,3)+F(I,4)*f(5,4) 

DO  9  1*1,  N 
DO  9  J-1,3 

9  F(l+9,J)”F(l,l)*F(6,j) 

DO  11  1=1,3 

11  F(l0,l43)  =  F(l,2)*F(6,l)4F(l,3)*F(7,l)+F(l,it)*F(8,l)4 

1F(1,5)*F(9,I) 

F( 10, 7 ) -F( 1,6) 

DO  13  1=1,3 
DO  13  J=l, 3 

13  F(I+10,  J+3)  =  F(I+1,2)*F(7, J)4F(l4l,3)*F(8,  J) +F( l4l,  4 ) * 
1F(9,  J) 

DO  15  1=1,3 

15  F(l+10,7)=F( 1+1,5) 
do  17  1=1,4 
DO  17  J-1,7 

17  F(I+13,  J)=F(I,1)*F(10,  J) 
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DO  19  1=1,  -4 

19  f(I+13j8)=f(I+9,7 ) 

DO  21  1=1,7 

21  F(l4, 1+8)=F(l,2)*F(lO,l)-+F(l,3)*F(ll,l)+F(l,4)*F(l2,l)  + 
1F(1,5)*F(13,I) 

DO  23  1=1,3 
DO  23  J=l,7 

23  F(l4l4,  J-+8)=F(l-tl,2)*F(ll,  J)4F(I+1,3)*F(12,  j)4F(l4l,4)* 
1F(13,  J) 

C  TRANSFER  THE  COEFFICIENTS  TO  B  MATRIX 
DO  40  1=1,  N 
DO  40  J=l,  15 

40  B(l,J)=F(l4l3,  J) 

C  ELIMINATE  VARIABLE  GAIN  CK3 
DO  42  1=1,3 
DO  42  J=l,  16 

42  B(l44, J)=B(l4l, J)*B(1,8) 

DO  44  1=2,4 
DO  44  J=l,  16 

44  B(l46, J)=B(1, J)*B(I,8) 

DO  46  1=1,3 
DO  46  J=l,  16 

46  B(ltlO, J)=B(lf4, j)-B(l+7, J) 

C  TRANSFER  COEFFICIENTS  TO  MATRIX  A  TO  SOLVE  EQUATIONS 
DO  48  1=1,3 
DO  48  J=l, 8 

48  A(I, J)=B(I410, J48) 

GO  TO  80 
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C  MANIPULATION  OF  E  MATRIX  TO  OBTAIN  COEFFICIENTS  OF  THIRD - 
C  ORDER  STATE  VARIABLE  EQUATIONS 

49  E(4,l)=E(l,5) 

E(5,1)=E(2,4) 

e(6,i)=e(3,4) 

DO  50  1=1,3 

50  E(l*6,l)=E(l,l)*E(4,l) 

DO  52  1=1,3 

52  E(I+6,2)=E(I+3,1) 

E(7,3)=E(1,2)*E(4,1)+E(1,3)*E(5,1)<-E(1,4)*E(6,1) 

DO  54  1=1,2 

54  E(I+7,3)=E(I+1,2)*E(5, l)*E(l+l,3)*E(6,l) 
do  56  1-1,3 
DO  56  J=l,3 

56  E(l+9, J)=e(i,1)*e(7, J) 

DO  58  1=1,3 

58  e(10,I+3)=e(1,2)*e(7,I)*e(1,3)*E(8,i)«-e(1,4)*e(9,i) 

DO  60  1=1,3 

60  E(l+9,7)=E(I+3,1) 

DO  62  1=1,2 
DO  62  J=l,  3 

62  E(I<-10,  J+3)=E(I+1,2)*E(8,  j)tE(lt-l,3)*E(9,  j) 

DO  66  1=1, N 
DO  66  j=l,7 

66  A(I,  J)=E(l4-9,  j) 

C  SOLUTION  OF  THIRD-ORDER  SIMULTANEOUS  EQUATIONS 

80  DO  82  1=1,2 
DO  82  J=l,  8 
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82  A(l+3, J)=A(I+1, j)*A(l,7) 

DO  84  1=2,3 
DO  84  J=l,  8 

84  A(H-4, J)=A(l, J)*A(I,7) 

DO  86  1=1,2 
DO  86  J=l, 8 

86  A(l+7, J)=A(l+3, J)-A(l+5, J) 

DO  87  J=l,  8 

87  A(10, J)=A(8, J)*A(9,5) 

DO  88  J=l,  8 

88  A(ll, J)=A(9, J)*A(8,5) 

DO  89  J=l,  8 

89  A(l2, J)=A( 10, j)-A(ll, J) 

C  SOLVE  FOR  VARIABLE  GAINS 

CKO=A (12, 8) /A (12, 6) 

CK1= ( A( 8, 8) -A (8, 6) *CKO)/(A (8, 4 ) *CKO-A (8,5)) 

CNK2=A (l,8)-A(l,4) *CK0*CK1-A (1,5) *CKI-A ( 1, 6) *CKO 

CDK2=A (l, 1 )*CK0*CK1+A ( 1 , 2) *CKlfA( 1,3) *CKO*A( 1, 7 ) 

CK2=CNK2/CDK2 

IF  (N.EQ.3)  GO  TO  91 

CK01=CK0*CK1 

CK02=CK0*CK2 

CK12=CK1*CK2 

CK012=CK0*CK1*CK2 

CKN3=B(1, 16) -B( 1,9) *CK012-B( 1,10) *CK12-B( 1,11) *CK02- 
1B( 1, 12 ) *CK01-B( 1, 13 ) *CK1-B( 1 , 14 ) *CKO-B( 1,15) *CK2 
CKD3-B(  1, 1 )  *CK012f  B(  1,2)  *CK12+B( 1, 3 )  *CK02+B(  1, 4 )  *CK01+- 
1B(  1,  5 )  *CK1  +  B(  1 , 6)  *CKO+B(  1,7 )  *CK2f  B(  1,8) 


(  •  )  (t  j:  r )  a  (  c.  A  A, 

(Y,l)A*(L,£)A-(T;,a+l)A  48 
s,£=i  58  oa 
8.X-T,  58  oa 

(t,,5+l)A-(T.,£4l)A=(t,,Y+I)A  58 

/ 

(5,e)A*(T,t8)'f -(V,OX)A  Y8 
(e,8)A*('Gie),\-(T.,£X)A  88 

,  '  .  8 ,Xrt  68  oa 

(t,£I>fcr(T,,0X)A=(it2l)A  68 
SHI  AO  3JQAIHAV  HO1?  3VIOS  0 
(5,21)  A\  (8 .  SX  )A*0X0 
((£,3)fl-OX''"'(t1,L  'll  )\{(  M5£'  CAO 

5,X)A-XHO*(5,X)A-Xa 
H0*(£,X)A+XX0*(2, X 


X6  OT  00  (E.PH.M)  ’ll 

£XO*£XO=2£XO 

2XO*XXO*OXO=2XO!IO 

-20X0* ( IX  t  X  )8-2IX3*( OX , X  )S-2X0X0* ( 6 , 1)3- ( 5l ,1)8 

t,x)a-oxo*<Ax,x)j  ■  o*(£x,x)a-xoxo*(2i,x)ax 

4  XOXO*  (  4  ,  X  )a+£OAO»  ( £  ,  1)3+21X0*  ( S  ,  X  )8+2£OXO*  ( X  ,1)8- £<E 

(8»x:  1  -  ‘  ;  ,.[)3+oxo*(5«x)a+ixo*(e,.r )ax 
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CK3  =  CKN3/CKD3 

WRITE  (6,90)  CKO,  CK1,  CK2,  CK3 

90  FORMAT  (26H  CONTROLLER  GAINS  CKO--CK3/4E16. 8) 

GO  TO  97 

91  WRITE  (6,92)  CKO,  CK1,  CK2 

92  FORMAT  (26H  CONTROLLER  GAINS  CKO--CK2/3E16 . 8) 

C  D(Z)  FOR  STEP  INPUTS 

C  ERRORS  FOR  THIRD- ORDER  SYSTEM 
ERO=A (l,8) 

ERl=A(l,8)-E(4,l)*CKO 

ER2=A(1,8)-E(7, 1)*CK0*CK1-E(7,2)*CK1-E(7,3)*CK0 
C  PLANT  INPUTS  FOR  THIRD- ORDER  SYSTEM 
HO=CKO*ERO 
H1=CK1*ER1 
H2=CK2*ER2 

WRITE  (6,94)  ERO, ER1, ER2 

94  FORMAT  (23H  ERROR  SIGNALS  ERO--ER2/3E16 . 8) 

WRITE  (6,95)  HO, HI, H2 

95  FORMAT  (26H  CONTROLLER  OUTPUTS  HO--H2/3EI6 . 8) 

CALL  EXIT 

C  ERRORS  FOR  FOURTH- ORDER  SYSTEM 
97  ERO=B(l,l6) 

ER1=B(1, 16)-F(5, l)*CKO 

ER2=B(1, l6)-F(6, 1)*CK0*CK1-F(6,2)*CK1-F(6,3)*CK0 
ER3=B(l,l6)-F(l0,l)*CK012-F(l0,2)*CK12-F(l0,3)*CK02- 
1F(10,4)*CK01-F(10,5)*CK1-F(10,6)*CK0-F(10,7)*CK2 
C  PLANT  INPUTS  FOR  FOURTH- ORDER  SYSTEM 


HOzCKO*ERO 
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H1=CK1*ER1 

H2=CK2*ER2 

H3=CK3*ER3 

WRITE  (6,98)  ER0,ER1,ER2,ER3 

98  FORMAT  (23H  ERROR  SIGNALS  ERO--ER3/4E16. 8) 

WRITE  (6,99)  H0,H1,H2,H3 

99  FORMAT  (26H  CONTROLLER  OUTPUTS  HO--H3/4E16 . 8) 

END 

$IBFTC  COEF 

SUBROUTINE  ARRAYS  (N,T) 

COMMON  A(12,8),B(13, 16) , E(l2, 7) , F(17, 15) 

C  THIS  SUBROUTINE  IS  USED  TO  PUNCH  THE  COEFFICIENTS  OF  THE 
C  STATE  VARIABLES  AND  THE  FINAL  VALUES  OF  THE  STATE  EQUATIONS 
C  INTO  THEIR  APPROPRIATE  MATRICES.  FOR  EXAMPLE,  THE  VALUES 
C  IMMEDIATELY  FOLLOWING  ARE  THE  REQUIRED  DATA  FOR  THE  PLANT 
C  G(S)=l/S(Sfl) (S+2) . 

A(l,8)  =  1.0 
A (2, 8) =0. 0 
A (3, 8)  =  0. 0 

E(l,l)=-(T/2.-3./4.fEXP(-T)-0.25*EXP(-2.*T) ) 

E(l, 2) =1. 

E(l,3)=0.5-0.5*EXP(-2.*) 

E(l,4)=0.5f0.5*EXP(-2.*T)-EXP(-T) 

E(l,5)=T/2.-0.75+EXP(-T)-0.25*EXP(-2.*T) 

E(2,l)=-0.5-0.5*EXP(-2.*T)fEXP(-T) 

E ( 2, 2 ) =EXP( -2 . *T) 

E ( 2, 3 ) =EXP( -T ) -EXP( -2 . *T ) 


■ 


'  lO 
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E(2,4)=0.5+0.5*EXP(-2.*T)-EXP(-T) 

E(3,1)=-1.+EXP(-T) 

E(3,2)=0.0 
E (3>  3 ) -EXP(-T) 

E(3,4)=1.-EXP(-T) 

RETURN 

END 


3 . 3  INPUT 


The  coefficients  of  the  state  variables,  in  terms  of 
the  sampling  period  T,  are  introduced  into  the  program  using 
the  subroutine  "ARRAYS'* .  This  allows  T  to  be  varied  and  the 
magnitudes  of  the  resulting  gains  to  be  examined.  Each 
coefficient  is  placed  in  a  specific  location  in  the  E  or  F 
matrix.  The  state  variable  equations  for  a  third-order 
system  with  the  coefficients  replaced  by  the  appropriate 
matrix  designations  are  shown  below. 


(k>  1  )T 


E(l,l)KkvE(l,2)  xx(kT)  t  E(l,3)x2(kT) 
i-  E(l,it)x3(kT)  t-  £(1,5)1^ 


(kpl)T 


=  E(2,l)Kkx1(kT)  P  E(2,2)x2(kT) 
P  E(2,3)x3(kT)  P  E(2,4)Kk 


(k+l)T 


=  E(3,l)Kkx1(kT)  P  E(3,2)x2(kT) 
+  E(3,3)x3(kT)  P  E(3,4)Kk 


(3-D 


(3-2) 

(3-3) 


The  coefficients  of  a  fourth-order  system  are 
similarly  read  into  the  F  matrix  locations. 


H*?  .  04?  .0  ;  ,S;3 


(  ••)'  -  ■  '  ’ 


(a1-  Hxs-.i*(ji,£)s 


TU3WI  £.£ 


:  .  "!  '  •'  .  .  .  _  If;  " 

gnieu  meigoiq  9riS  oini  fc9ou&oT:;fnJ:  arts  .T  bol/i9q  80iXqni88  srij 

;a  erii 
13  esbu^lngsffl 

to  H  srfct  ni  noll&ooL  oillosqa  e  nl  b9osIq  si  cfnsi oillsoo 
ci9b‘X0-5*xirt;t  6  iorl  an 

'  •  “  "  '  ’  .  :  '*  r  irf  r\-  :  ■  ■■  i 


-  -■  .  •  I  .  ■  '■  -  *  1 


>  i(c,i)2  ■*  ’>i  :;  •<[;  ~i  t 


. 

' 1  u 


.  r  ►'.£'  *  (  '■))•<(-,  4 


•  ■  ' 


,  .  ,  rr-t  j  .  .  .  .  . 


28  - 


The  predetermined  final  values,  to  which  the  state 
equations  are  equated,  are  placed  in  the  A  or  B  matrix  as 
follows : 

for  a  third-order  system 

A(l,8)=  value  of  x^  at  time  3T 

A(2,8)=  value  of  Xg  at  time  3T 

A(3j8)=  value  of  x^  at  time  3T 

for  a  fourth-order  system 

B(l,l6)=  value  of  x-^  at  time  4T 

B(2,16)-  value  of  Xg  at  time  4t 

B(3*l6)=  value  of  x^  at  time  4t 

B(4,l6)-  value  of  x^  at  time  4t 

To  complete  the  input  information,  the  numerical 
values  of  the  order  (N)  of  the  system  and  the  sampling  period 
(T  seconds)  are  punched  on  data  cards.  The  last  data  card 
has  a  9  punched  into  the  location  for  system  order.  When  the 
computer  reads  this  card,  the  execution  of  the  program  is 
stopped . 

3.-4  OUTPUT 

Shown  below  is  a  sample  print-out  of  the  program  for 
the  plant  with  transfer  function  G(s)=l/s(s+l) (s+2) . 

3  0.100E  01 

CONTROLLER  GAINS  CK0--CK2 

0.36591675E  01  -0.26591258E  01  0.26386268E  01 

ERROR  SIGNALS  ER0--ER2 

0 . 10000000E  01  0.69246299E  00  0.69042742E-01 


[ 
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CONTROLLER  OUTPUTS  H0--H2 
0.36591675E  01  -0. 18413462E  01  0.18217802E  00 

The  first  line  of  the  output  is  a  print-out  of  the 
input  data  card.  The  integer  3  is  the  order  of  the  system 
and  the  exponential  number  0.100E  01  is  the  sampling  period 
in  seconds,  one  second  in  this  case. 

The  CONTROLLER  GAINS  CK0--CK2  are  the  variable  gains 
Kq  =  3.6591675 
K-l  -  -2.6591257 
k2  =  2.6386267 

The  ERROR  SIGNALS  ER0--ER2  are  the  errors  at  the 
sampling  instants 

e(0*)  =  1.00000000 
e(T*)  =  0.69246299 
e ( 2Tf ) =  0.069042742 

The  CONTROLLER  OUTPUTS  H0--H2  are  the  inputs  to  the 

plant 

h(0^ )  =  3.6591675 
h(T+)  =-1.8413462 
h ( 2T1' )  -  0.1817802 
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CHAPTER  IV 


TRANSFER  FUNCTIONS  OF  DEADBEAT  CONTROLLERS 

Using  the  program  described  in  the  previous  chapter, 
the  transfer  functions  of  deadbeat  controllers  for  various 
systems  will  be  obtained  and  the  results  confirmed  on  the’ PACE 
231-R  analog  computer.  The  procedure  used  to  solve  a 
system  containing  a  zero  in  the  transfer  function  will  also 
be  described. 

4.1  THIRD- ORDER  SYSTEM  WITH  SIMPLE  POLES 

For  a  plant,  which  has  a  transfer  function  consisting 
of  simple  poles  only,  the  deadbeat  response  that  is  obtained 
has  no  overshoot  irrespective  of  the  sampling  period.  But, 
the  sampling  period  does  affect  the  magnitudes  of  the 
variable  gains.  Considering  the  plant  G(s)=  l/s( s+l) (s+2) 
and  its  flow  graph  shown  in  Figure  2-6,  solving  for  sampling 
periods  T=1  second  and  T=0.1  second,  yields  the  transforms 

T=1  second  n( ,  )=3- 6591675-1 .  fltlB^z-i+O.  l8217802z'2 

1 . 0000000+0 . 692i|6299z_1+0 . 0690^t27itz‘2 

T=0.1  second  n( *)=1159 •  4171-1998. 33^z-l+858.917l8z-2 

1. 0000000+0. 820 sgilSSz-i+O.  I5it4l565z'2 

The  validity  of  the  values  for  T=1  second  was  verified 
on  the  analog  computer  using  the  rep-op  method*.  The  analog 


A  description  of  this  method,  developed  by  Professor 
Y.  J.  Kingma,  may  be  found  in  Appendix  A  of  the  master’s  thesis 
Sushil  K.  Sarna  on  ' Sampled-Da ta  Conditional  Feedback  Systems', 
University  of  Alberta,  September,  1966, 
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FIGURE  4-1.  ANALOG  SIMULATION  FLOW  GRAPH  FOR  G(s)=_, - =4, - - 

s(s+l) (s+2) 


FIGURE  -4-2.  STATE  VARIABLES  OF  COMPENSATED  SYSTEM 


, 
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diagram  is  shown  in  Figure  4-1  and  the  results  in  Figure  4-2. 
The  analog  diagram  is  time-scaled  for  optimum  display  on  a 
multi-channel  oscilloscope. 

Figure  4-2  shows  that  the  output  (A)  has  reached  a 
final  value  after  three  sampling  periods  with  no  overshoot 
and  the  other  state  variables  (B  and  C)  have  become  zero  at 
the  same  instant. 

As  the  duration  of  the  sampling  period  is  decreased, 
the  magnitudes  of  the  variable  gains  and  the  plant  inputs 
are  increased.  This  is  expected,  as  the  output  must  rise  to 
its  final  value  in  a  shorter  time.  The  z-transform  of  the 
digital  controller,  for  a  sampling  period  of  T=0.1  seconds, 
exhibits  the  large  variable  gains  and  plant  inputs  required 
as  compared  to  those  for  a  sampling  period  of  T=1.0  seconds. 
Thus  in  practice,  the  gains  that  can  be  realized  by  control 
equipment  may  determine  the  sampling  period. 

4.2  EFFECT  OF  REARRANGING  SYSTEM  FLOW  GRAPHS 


Using  an  example,  the  system  flow  graph  of  Figure  2-6 
is  rearranged  to  that  of  Figure  4-3. 
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The  state  equations  written  from  the  flow  graph  of 
Figure  4-3  are 


xl(s)  -  -Kks  ^ ^+s  1  x1(tQ)  +  s"2(lf2s"1)/^ 


MV 


Xo(s)  = 


X3(s) 


3  +  r(t0)Kk  C3"4/^] 

-Kks'3(lfs‘1)/^jx1(t0)  <■  [s‘-jx2(t0) 

+  S‘2(l4-s'1)/^jx3(t0)  ¥  r(t0)Kk  js*3(i*s-l)/^J 

r  j2Kk8’2(lts’1)'&  xl(to)  +  [s-1(l^s-1)^jx3(t0) 

*■  s-2(l*s'1)/^  rCtg)!^ 


(4-1) 


A  1  ( 1+3S--U2S  2) 


Comparing  equations  4-1  with  equations  2-4,  it  is 
noted  that  the  coefficients  of  the  state  variables  are 
different  although  the  plant  has  the  same  transfer  function 
To  solve  equations  4-1,  the  final  values  of  the  state  vari¬ 
ables  also  differ  from  those  used  to  solve  equations  2-4. 
The  final  values  for  both  sets  of  state  variables  are 
Figure  2-4  Figure  4-3 

Xl(3T)  =  1.0  x1(3T)  =  1.0 

x2(3T)  =  0.0  x2(3T)  =  1.0 

x3(3T)  =  0.0  x3(3T)  =  0.0 


The  transfer  functions  of  the  digital  controller, 
obtained  using  the  digital  program  and  a  sampling  period  of 
T=1.0  seconds,  are  identical  in  both  cases.  This  shows  that 
the  method  of  drawing  the  flow  graph  does  not  affect  the 
results  obtained  for  the  digital  controller. 
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4.3  FOURTH- ORDER  SYSTEM  WITH  SIMPLE  POLES 

For  the  plant  with  transfer  function  — , - 1 _ r_ _ _ 

s ( s+1 ) ( s+2 ) ( s+3 ) 

and  a  sampling  period  of  T=1  second,  the  z-transform  of  the 
digital  controller  was  found  to  be 

D(z)  =  13.553  -  6.3886z~^  4-  0.86460z~|~  -  Q.028636z~^ 
1.0000  +  0.83846z_1f  0.19953z"2  t  0.0044535z"3 


With  the  aid  of  the  analog  simulation  flow  graph  of 
Figure  4-4,  the  transfer  function  was  tested  by  the  rep-op 
method  and  the  results  are  shown  in  Figure  4-5. 


FIGURE  4-4.  ANALOG  SIMULATION  FLOW  GRAPH 


4.4  THIRD-ORDER  OSCILLATORY  SYSTEM 


The  transfer  function  of  a  third-order  oscillatory 

plant  is  of  the  form  ,  p . .  The  analog  simulation 

flow  graph  of  the  plant  for  a  damping  ratio  of  ^  a 
natural  frequency  of  Ljn=  TT ,  and  a  period  of  oscillation  =  2 
seconds  is  shown  in  Figure  4-6.  Transforms  of  the  digital 
controllers  were  obtained  using  the  digital  program  and  tested 
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FIGURE  4-5.  STATE  VARIABLES  OF  FOURTH-ORDER  SYSTEM 


using  the  rep-op  method  for  sampling  periods  of  T=1  second 
and  T=2  seconds.  The  results  are  shown  in  Figures  4-7  and 
4-8  for  the  following  z-transforms 

T=1  second:  D(z)  -  6.9373  +  2.6325Z'1  +  0.29979z~2 

1.0000  +  0. 57^272' x+  0.071523z^ 

T-2  seconds:  D(z)  -  5-2259  -  0. 30087*- j  +  0.0097588z-2 

1.0000  +  0.10156z_i  +  0. 0005769IZ"2 


FIGURE  4-6.  ANALOG  SIMULATION  FLOW  GRAPH  FOR  G(s)=  — — r^T 

S  S  4-  £_ Jj  ( As  S  / 

Choosing  a  sampling  period  close  to  the  time  period 
of  the  oscillating  system,  results  in  slight  overshoots  as 
can  be  seen  in  Figure  4-8,  but  the  output  has  almost  reached 
the  value  of  the  input  in  one  sampling  period.  The  remaining 
two  sampling  periods  are  required  to  level  off  the  output  to 


•.  :  O'*  3 :■  t*  q  3x1.: .  ..  \  c  nr  <rU 


ri  -ori  v  +n'i  •  if  ilt/r  >**  tme^^\r8  CJtnao  o 

. 


37 


FIGURE  *1-7 . 


STATE  VARIABLES 


OF  OSCILLATORY  SYSTEM  FOR 


T=  1  SECOND 


FIGURE  4-8.  STATE  VARIABLES  OF  OSCILLATORY  SYSTEM  FOR 


T=  2  SECONDS 
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the  designated  value.  For  a  sampling  period  of  one-half  the 
time  period,  the  output  is  more  stable  as  there  is  no 
tendency  to  oscillate  about  the  final  value.  The  value  of 
the  output  in  Figure  4-7  after  two  sampling  periods  is 
approximately  equal  to  the  value  of  the  output  in  Figure  4-8 
after  one  sampling  period,  but  the  time  interval  is  the  same 

A 

in  both  cases.  Thus,  a  knowledge  of  the  period  of  oscillation 
is  helpful  in  choosing  a  sampling  rate  that  will  control  the 
output.  For  higher  order  systems  with  unfactored  denominators, 
a  program  for  solving  the  roots  of  a  polynomial  is  described 
in  Chapter  V.  Once  the  roots  have  been  found,  the  period  of 
oscillation  is  easily  determined. 

4.5  THIRD- ORDER  SYSTEM  WITH  A  ZERO 

A  third-order  plant  containing  a  zero  has  a  transfer 

function  of  the  form  G(s)  = - r  .  When  drawing  the 

s(s+b)(s+c) 

flow  graph  for  this  system,  an  extra  node  is  required  to 
add  the  effect  of  the  zero  onto  the  output.  The  open-loop 
flow  graph  may  be  drawn  as  in  Figure  4-9a  or  in  Figure  4-9b. 


-c 


a.  Method  1 


b.  Method  2 


FIGURE  4-9.  OPEN- LOOP  FLOW  GRAPHS 


Applying  Mason's  gain  formula  to  both  flow  graphs 
results  in  the  same  transfer  function,  but  Method  1  cannot  be 
set  up  directly  on  the  analog  computer  as  an  additional  ampli- 
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fier  is  required  to  add  the  ',aX2n  term  to  the  "X1'f  term. 
Therefore,  it  is  the  value  of  Y  and  not  X2  that  is  subtracted 
from  the  input  at  the  sampling  instants.  The  complete  flow 
graph,  including  a  digital  controller  is  shown  in  Figure  4-10. 


To  solve  the  system  of  Figure  4-10  requires  four 

equations,  one  for  each  node,  but  the  equation  for  Y  is  a 

linear  combination  of  X^  and  Xg,  thus  there  are  three 

independent  equations  in  three  unknowns. 

The  program  of  Chapter  III  cannot  be  used  to  solve 

these  equations  because  of  the  added  equation  for  Y.  Solving 

the  equations  for  the  plant  G(s)=  — , ......  and  T=1  second, 

s  ( s  + 1 )  (  s  +  2 ) 

results  in  the  D(z)  =  3 »  6392  - — 1_:  84l3z - +  0  •  l82l8z — 

1.0000  +  0. 32693z_1-  0.065428z"2 


which  has  the  same  numerator  as  the  D(z)  obtained  for 

G(s )-  ^  ,  that  is 

v  '  s(s+l) (s+2) 
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D(z)  =  3.6592  -  1.84l3z  1  4-  0.  i82l8z  2 

1.0000  +  0.692-46Z"1  +  0.c69043z"'d 

Thus  if  the  zero  term  were  neglected,  the  program  of 

Chapter  III  could  be  used  to  solve  for  the  inputs  to  the 

plant,  h(0+),  h(T+)  and  h(2T+).  Applying  these  inputs  to  the 

open-loop  plant  and  measuring  the  difference  between  the 

final  output,  magnitude  assumed  as  one,  and  the  output  at  the 

sampling  instants,  the  error  signals,  e(T+)  and  e(2T+‘),  can 

be  calculated. 

Using  the  rep-op  method  with  the  flow  graph  of 

Figure  4-11  and  the  above  plant  inputs,  h(nT+),  for 

G(s)=  1  ,  the  outputs  for  G-i(s)  =  1  , 

s(s+l) (s  +  2)  s(s+l) (s*2) 

Go(s)=  S4-.5  and  Gq(s)=  s+1.5  were  observed, 
s ( s+l) (s+2)  s(s+l) ( s+2) 

The  graphs  of  the  outputs  are  shown  in  Figure  4-12  for  a 

sampling  period  of  one  second. 


FIGURE  4-11.  ANALOG  FLOW  GRAPH  FOR  SYSTEMS  WITH  A  ZERO 

All  three  outputs  have  become  deadbeat  after  three 
sampling  periods,  but  the  further  the  zero  is  removed  from 


to  ms'igoTq  erl$  tbe3o9f$su  rmsj  oiss  9rfcf  *tl  aurfT 

< 

sricJ  ol  eJuqnl  9riJ  not  sv  os  otf  oeeu  r*d  b£uoo  III  rfsctQsriO 

srict  3&  JuqJuo  oriJ  bns  t9f}0  as  bs^L/sas  abucElngs/n  t;fuq:fuo  IenJt'1 

.  E  9  3  Bltro  l  BO  B<1 

. 

.bnoosa  .-.'no  ‘  o  coiisq  gnllqmee 


OHSS  \  HTIVf  3M3T3Y8  HC  I  H3AJ  0  W0J3  DO  A  ..  -H  g-TODM 

mori  bavonm  ajt  oiss  erfi  isrlchtt/'  en’^  j./ci  taboi  teq  Ilqmss 


the  origin  in  the  s-plane,  the  larger  the  overshoot  becomes. 
For  systems  with  the  pole  and  zero  separated  by  a  decade 
or  more,  no  overshoot  will  occur  in  the  output. 


FIGURE  4-12.  OUTPUTS  FOR  SYSTEMS  WITH  A  ZERO  COMPARED  TO 


SYSTEM  WITH  NO  ZERO 


CHAPTER  V 


COMPUTATION  OF  ROOTS  OF  A  POLYNOMIAL  USING  ROUTH’ S  STABILITY 
CRITERION 


The  concept  and  method  of  applying  Routin' s  stability 
criterion  to  computing  the  roots  of  a  polynomial  was  first 
introduced  by  Professor  Y.  J.  Kingma  and  a  subsequent  program 
was  written  by  Professor  D.  H.  Kelly  of  the  Electrical 
Engineering  Department  at  the  University  of  Alberta.  The 
object  of  this  thesis  was  to  rewrite  the  original  program, 
written  in  FORTRAN  II  computer  language,  into  FORTRAN  IV 
language,  double  precision  the  new  program  and  check  the 
accuracy  obtained  using  the  program. 

5.1  SCOPE  OF  DIGITAL  PROGRAM 

The  Routin' s  stability  criterion  is  used  in  control 
system  analysis  to  determine  the  number  of  roots  of  the 
characteristic  equation  of  a  closed-loop  system  that  have 
positive  real  parts.  This  is  accomplished  by  computing  the 
Routh  array  and  counting  the  number  of  sign  changes  in  the 
first  column  of  the  array.  Each  sign  change  corresponds  to 
a  root  with  a  positive  real  part.  A  description  of  Routh' s 
stability  criterion  may  be  found  in  the  literature^. 

The  program  finds  the  roots  of  a  polynomial  by 
successively  shifting  the  imaginary  axis  and  calculating  the 
Routh  array.  The  axis  is  shifted  by  increasingly  smaller 
increments  until  it  converges  on  the  real  part  of  the  root 
to  the  accuracy  specified  in  the  input  data.  Roots  of  any 
order  of  multiplicity  may  be  located  using  the  program. 
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If  there  are  no  more  than  two  complex  pairs  with  the  same 
real  part,  the  imaginary  parts  of  the  root  are  then  calcu¬ 
lated.  If  there  are  three  or  more  pairs  of  complex  roots  with 
the  same  real  part,  the  coefficients  of  an  auxilliary  equa¬ 
tion  are  printed  out  which  may  be  inserted  back  into  the 
program  as  new  data  to  solve  for  the  imaginary  parts. 

5.2  LIST  OF  TERMS  USED  IN  PROGRAM 


N  -  degree  of  polynomial 

A(j)  ---  coefficients  of  polynomial 

AA(j)  --  coefficients  of  polynomial  with  shifted  axis 

2  N 

C(l,J)  -  array  of  coefficients  of  (s  +  l),  (s+1)  ,  . (sfl) 

B(l,J)  -  coefficients  of  Routh  array 

NR  -  number  of  roots  still  to  be  calculated 

AL  -  distance  axis  is  shifted  from  original  position 

NRR  -  number  of  roots  with  a  real  part  located  at  AL 

ALO  -  size  of  initial  rightward  step 

ANO  -  size  of  initial  leftward  step  (negative) 

El  -  value  used  to  test  B(l,l)  for  zero 

E2  -  value  that  replaces  B(l,l)  if  B(l, 1)  <  El 

E3  -  desired  accuracy,  equal  to  value  of  Incremental  shift 

value  of  real  part  of  root 


5.3  CALCULATIONS  PERFORMED  USING  COMPUTER  PROGRAM 

This  program  is  dimensioned  to  accomodate  polynomials 
of  degree  sixteen.  For  higher  order  polynomials,  modify  the 
dimension  statement  to 

DIMENSION  B(N2,N2),  A(N1),  AA(N2),  C(N2,N2) 
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where  N2=order  of  polynomial  +  2 
Nl=order  of  polynomial  4-  1 

Given  a  polynomial  of  the  form 

A(l)sN  +  A(2)sN_1  4  .  +  A(N4l)  . (5-1) 

the  first  step  is  to  shift  the  imaginary  axis  of  the  s-plane 
to  the  left  by  an  amount  ,  ANO,  and  calculate  a  new  poly¬ 
nomial  of  the  form 

A ( 1)  ( Sfo()N  ^  A(2)(s+o<)N"1  +  .  +  A ( N+l )  . (5-2) 

which  reduces  to 

AA(1)sN  +  AA(2)sN"1  4  .  +  AA(N+1)  . (5-3) 

The  reduction  of  equation  (5-2)  to  equation  (5-3)  is 
carried  out  with  the  aid  of  an  array,  C(l,j),  containing 

the  coefficients  of  (s+1),  (s+l)^,  . .  (s+l)^  as  shown  in 

Figure  5-1 • 
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FIGURE  5-1.  C- ARRAY 

The  first  row  contains  the  coefficients  of  (s+-l),  the 

2 

second  row  the  coefficients  of  (sfl)  ,  and  so  on  down  to  the 
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N-th  row  which  contains  the  coefficients  of 


(s*l) 


N 


The  coefficients  of  equation  (5-3)  are  formed  thusly: 
for  N  =  2,  AA ( 1 )  =  A(l) 

AA(2)  =  A(l)C(2,2)*  +■  A(2) 

AA(3)  =  A(l)C(2,3)o<2  *  C(l,2)A(2)<v  -  A(3) 


for  N  =  3, 

AA  ( 1 )  =  A  ( 1) 

AA(2)  =  A(l)C(3,2)of  *  A(2) 

AA(3)  =  A(l)C(3,3)o(2*  A(2)C(2,2)ot  +  A(3) 

AA (4 )  =  A(1)C(3,4)o(3+  A (2 ) C (2, 3 )o(2+  A(3)C(l,2)o{  t  A(4) 
and  for  N  =  4, 


AA  ( 1 )  =  A  ( 1 ) 

AA ( 2)  =  A(1)C(4,2)o<  -  A (2 ) 

AA  ( 3 )  =  A(l)C(4,3)»<2tA(2)C(3,2)0<  +  A(3) 

AA(4)  =  A(l)C(4,4)0<3+  A(2)C(3,3)o(2+  A(3)C(2,2)s<  +  A(4) 

AA ( 5 )  =  A(l)c(4,  5)o/*  +  A(2)c(3,4)^3<-  A(3)C(2,3K2  + 

♦  A(4)C(l,2)o{  +  A(5) 

For  higher  values  of  N  the  coefficients  are 
calculated  in  a  similar  manner. 

The  coefficients  of  equation  (5-3)  form  the  first  two 
rows  of  the  Routh  array,  B(l,J),  as  shown  in  Figure  5-2.  The 
remainder  of  the  Routh  array  is  built  up  using  equation  (5-^). 
B( I, J)  =  B( 1-2, J+1 )  -  B(I-1, J+l)*B(l-2,l)/B(l-l,l). . . (5-4) 

After  each  row  is  formed,  the  first  element,  B(l,l), 
is  tested  for  a  ’’zero"  value.  That  is,  is  B(l,l)<El,  where 
El  is  a  very  small  number?  If  it  is  less  than  El,  each 
remaining  element  in  the  row  is  tested  for  zeros,  and  if  any 
of  these  have  a  value,  B(l,l)  is  replaced  by  a  small  number  E2. 


(S) A  4  *(2,2)0(X)A  -  (2)AA 
(E)A  .4  >>(2)A(2,X)0  4  S»(£,S)0(X)A  =  (£)AA 


(S) A  4  y5(S,£)0(X)A  =  (S)AA 
(e)A  4  *(2,2)0(2)A  ♦*(£, £)0(X)A  =  {  )AA 


(|i)A  *  jo(S  tI  )0( 6) A  ^S>o(e,S)0(S)A  >  AA 


(2) A  4  )o(;a,t1)0(X)A  *  (2) AA 
(£)A  4  s.(2,£)0(S)A4>,(£,t‘)0(X)A  *  (£)AA 

4  %(£,2)0{£)A  4^«1,£)0{2)A  4  a(5,H)0(i)A  «  (5)AA 
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5 

FIGURE  5-2.  ROUTH  ARRAY  [b(I,J)] 

If  the  entire  row  is  zero,  the  coefficients  are  replaced  by 
the  coefficients  of  the  derivative  of  the  previous  row. 

Continuing  this  process  of  forming  a  row,  testing 
it  and  making  any  necessary  changes  the  complete  Routh  array 
is  constructed.  The  sign  changes  in  the  first  column  are 
counted  and  thus  the  number  of  roots  to  the  right  of  the 
shifted  axis  is  determined.  If  this  is  equal  to  the  order 
of  the  polynomial,  a  search  is  started,  otherwise  the  axis 
is  again  shifted  to  the  left  and  the  procedure  repeated 
until  all  roots  are  to  the  right  of  the  axis. 

The  search  is  started  by  shifting  the  axis  to  the 
right  by  an  amount  ALO,  and  the  calculations  are  carried  out 
as  above.  If  the  number  of  roots  to  the  right  has  not 
changed,  the  axis  is  again  shifted  to  the  right,  otherwise 
an  incremental  shift  to  the  left  is  made.  Right  and  left 
shifts  of  decreasing  size  are  made  until  the  axis  converges 
on  the  real  part  of  the  root  to  the  desired  accuracy  E3. 

The  number  of  roots  and  AL,  the  real  part  of  the 
root,  are  printed  out.  The  number  of  roots  found  is  subtracted 
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from  the  order  of  the  polynomial  and  the  result  stored  in 
a  location  NO.  This  is  the  value  that  is  now  used  for  a 
comparison  when  counting  the  number  of  sign  changes  in  the 
first  column  of  the  Routh  array.  If  the  number  of  sign 
changes  is  equal  to  NO,  a  search  is  again  started.  The 
search  is  carried  on  until  all  the  roots  are  found,  that 
is  when  NO  =  0. 

If  there  are  2  to  5  roots  with  the  same  real  part, 
there  can  be  2  or  4  complex  roots  which  the  program  will 
calculate  from  the  coefficients  of  the  L-th  row  of  the 
Routh  array  where 

L  =  order  of  equation  -  number  of  roots 
If  there  are  6  or  more  roots  with  the  same  real  part,  the 
coefficients  of  an  equation  whose  roots  are  the  imaginary 
components  desired  are  printed  out.  The  equation  may  be 
solved  using  the  program  and  the  coefficients  as  input  data. 

5 . 4  INPUT 

The  first  data  card  contains  the  values  of  N,  El,  E2, 
E3,  ALO  and  ANO  in  their  respective  order.  The  field 
specification  for  punching  the  numerical  values  on  the  card 
is  FORMAT  (I5,5D10.2). 

A  choice  of  0.01  for  E2  can  be  used  when  solving 
any  particular  polynomial  as  this  value  is  not  critical. 

The  values  for  E3,  the  desired  accuracy,  and  El,  the  value 
used  in  the  zero  test,  are  related  in  magnitude,  E3  being 
the  smaller.  Although  E3  is  the  desired  accuracy,  which  when 
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achieved  terminates  the  calculations,  both  El  and  E3  are  used 
in  tests  designed  to  terminate  the  iterative  calculations. 

If  DA/AL  <  E3,  where  DA  is  the  magnitude  of  the  incremental 
step  and  AL  the  distance  the  root  is  from  the  original 
imaginary  axis,  then  the  specified  accuracy  has  been  obtained 
and  the  real  part  of  the  root,  AL,  is  printed  out.  This  is 
the  dominant  test  for  large  roots,  but  if  the  root  is  very 
close  to  the  original  axis,  that  is,  AL  very  small,  then 
DA  <  El  is  the  dominant  test.  This  test  ensures  that  roots 
less  than  1  are  correct  to  a  specified  number  of  decimal 
places . 

—  HO  ft 

For  example,  choosing  E3  =  10”  and  El  =  10~  will 

ensure  that  roots  with  real  parts  less  than  one  will  be 

correct  to  eight  places  after  the  decimal.  Also,  roots 

with  real  parts  less  than  100  and  larger  than  one  will  be 

correct  to  eight  places  after  the  decimal  and  to  an  accuracy 
-8 

of  10~  per  cent. 

The  values  selected  for  ANO  and  ALO,  the  initial 
leftward  and  rightward  axis  shifts,  have  a  bearing  on  the 
computation  time.  The  choice  for  ANO  should  be  approximately 
three  times  as  large  as  that  for  ALO.  As  an  example,  for 
roots  at  s  =  -900,  -1,  ANO  =  -3,  and  ALO  =  1,  300  iterations 
are  required  to  reach  the  root  at  -900  and  another  900 
iterations  to  reach  the  root  at  -1.  A  choice  of  ANO  =  -30 
and  ALO  =  10  reduces  the  iterations  to  30  and  90  respectively, 
thus  speeding  up  the  solution.  The  saving  in  time  is  only 
a  few  seconds  if  the  program  is  being  run  on  the  IBM  7040 
computer,  but  is  several  minutes  when  using  the  IBM  1620 


«cis  Co  -  !  I.  <00*?  ed  t  5.  ■: ':b  art*.  !•  >  ■ 


©d  rn»  sno  nsr^  aa el  a^isq  lee* ^  w  a;*00'  £ 

3jo©*x  ,oe£A  .Ismioeb  9ri$  istfla  eeoBlq  Jrigls  otf  ctoi  'lo-: 

■ 


,01  ,-3C«n»9  ns  «A  .C  A  ':  >1  3-sd3  u  -«ib:  a  »  13 


,  sc  -13  .  T  ..  Jo  1  «  i  -> 


OP-  =  o  it  °fO  £>o  C”  A  .!-'•■  o  f  *r  iioBe><  oJ  &no  i'  oj 


.3  oi  -ns  C:  o  c-  .'  vre;  t  srtt  ieoulie*  01  0.!  be 


%i  io  a l  eml*  nl  ai  :  vst  ertT  .noJUuioa  erlcT  qu  gnibesqa  aurtt 


-  49  - 


computer  which  is  60  times  slower. 

By  studying  the  coefficients  A(2)  and  A(Nfl)  of  the 
given  polynomial,  the  user  of  the  program  may  obtain  an 
approximation  to  the  order  of  magnitude  of  the  roots.  The 
coefficient  A(2)  is  the  sum  of  the  roots  of  the  polynomial 
and  A(N+l)  is  the  product  of  the  roots. 

The  remaining  data  cards  contain  the  coefficients  of 
the  polynomial.  They  are  punched  onto  the  cards  in  the  order 
A(l),  A(2),  ....,  A(N+l)  according  to  the  field  specification 
FORMAT  (4Dl8. 10) . 

5.5  FORTRAN  IV  COMPUTER  PROGRAM  FOR  ROOT  SOLUTION 

The  flow  chart  for  this  program  is  in  Appendix  B. 

The  program  is  interspersed  with  COMMENT  statements  which 
indicate  the  calculations  that  are  to  follow. 


$IBFTC  RCM 

DOUBLE  PRECISION  B,A,AA,C,E1,E2,E3,AL0,AN0,DA,H,AL,R1, 
1R2,R3,R4 

DIMENSION  B( 18, 18) , A( 17) , AA( 18) , C(l8, 18) 

1  READ  (5,80)  N,E1, E2, E3,  ALO, ANO 
80  FORMAT  (15,5^10.2) 

Nl=N-*-l 

N2=N+2 

AL=ANO 

KM1=0 

NN=0 


NR=0 


J  1  ■  ■  :"»*> 
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KK=0 

KS=0 

DA=AL0/2.D0 
C  ZERO  OUT  AA  ARRAY 
DO  2  1=1, N2 

2  AA (l)=O.DO 

C  CALCULATE  COEFFICIENTS  OF  (S+l),  (S+l)**2,  . . . . (S+l)**N 
C ( 1 , 1 ) = 1 . DO 
C ( 1 , 2 ) = 1 . DO 
DO  3  1=2, N 
C(l, l)=l.DO 
C(l,Ifl)=l.DO 
DO  3  J=2, 1 

3  C(l,j)=C(l-l,j)-C(l-l,J-l) 

PE=E3*100 . DO 

WRITE  (6,81)  N, PE 

81  FORMAT (22H1  ORDER  OF  POLYNOMIALS I3/16H  PERCENT  ERROR 
ID  10. 2/) 

READ  (5, 82)  (A ( I) , 1=1, N1 ) 

82  FORMAT  (4Dl8. 10) 

WRITE  (6,83)  (A ( i) ,  1=1, Nl) 

83  FORMAT  (1X,4D17.10) 

C  ZERO  OUT  ROUTH  ARRAY 

AA(1)=A(1) 

DO  5  1=1, Nl, 2 
JJ=N2/2-l/2fl 
DO  5  J=1,JJ 
B( I, J) =0 . DO 


oa.x=(x,x)o 

Lfoy  oa.c*(a,x)o 


fij  .  :  ■■  •  rr 


oa.x=(£tj)o 

(l-t,X-l)0-(X.,I-X)0=(Lvl)0  e 

o  m  TM20H-R  HdX\£I  ,=J.\3.'.0:iYJ0‘;  20  flaano  XHSS)T  =**«  & 

(OX.Sraii)  TAMHOH  s8 
;  u  .-xa«i<xx)  'crswo'i  £8 


5  B(I+1, J)=O.DO 
GO  TO  44 
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C  FILL  IN  FIRST  2  ROWS  OF  ROUTH  ARRAY 
40  NN-0 
NR=0 

DO  6  1=1, Nl, 2 
J=(l+2)/2 

b( 1, j)=aa( i) 

6  B(2, J)=AA(I+1) 

C  TEST  SECOND  ROW  FOR  ZERO  COEFFICIENTS 
K=2 

JJ=N2/2 

KM1=K-1 

NN=N2-K 

IF(DABS(B(K, 1 ) ) .GE.E1)  GO  TO  19 

10  DO  11  J=2, JJ 

IF(DABS(B(K,J)).GE.E1)  GO  TO  13 

11  CONTINTJE 

DO  12  J=1,JJ 
H=NN+2-2*J 

12  B(K, J)=B(KM1, J)*H 
GO  TO  19 

13  B(  K,  1)  =E2 

19  IF(B(K,1)*B(KM1,1) .GE.O. )  GO  TO  15 

14  NR=NR+-1 

C  FILL  IN  NEXT  ROW 

15  DO  7  1=3, Nl, 2 


(X)AA*(T>  ,1  )3 
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JJ=N2/2-l/2 
DO  8  J=1,JJ 

B(I,  j)=B(l-2,  J<-l)-B(l-2, 1)*B(I-1,  Jfl)/B(l-1, 1) 

8  CONTINUE 

C  TEST  FOR  ZERO  COEFFICIENTS 
K=I 

KM1=K-1 

NNrN2-K 

IF(DABS(B(K,I)).GE.E1)  GO  TO  79 

70  DO  71  J=2,JJ 

IF(DABS(B(K, J) ) .GE.El)  GO  TO  73 

71  CONTINUE 

DO  72  J=1,JJ 
H=NNf 2-2* J 

72  B(K,J)=B(KM1,J)*H 
GO  TO  79 

73  B(Kj l)=E2 

C  TEST  FOR  SIGN  CHANGE 

79  IF(B(K,l)*B(KMI,l).GE.O.)  GO  TO  75 

74  NR=NR+1 

C  FILL  IN  NEXT  ROW 

75  DO  9  J=1,JJ 

B( If  1, j)=B(l-l, J+l)-B(l-l, 1)*B(I, Jf 1)/B(1, 1) 

9  CONTINUE 

C  TEST  FOR  ZERO  COEFFICIENTS 
K=I+1 


KM1=K-1 


1~y  »T0I*.  -.300  Of.  H03  1  i3T 

1-3 
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£Y  ot  oo  ( .o.30.(x,i!<ra)a*(x,H)a)?i  eY 


53  - 


NN=N2~K 

IF(DABS(B(K, l)).GE.El)  GO  TO  99 

90  DO  91  J=2, JJ 

IF(DABS(B(K, J) ) .GE.El)  GO  TO  93 

91  CONTINUE 

DO  92  J=1,JJ 
H=NN+2-2*J 

92  B( K, J) =B( KM1 ,  J ) *H 
GO  TO  99 

93  B(K,  1)=E2 

C  TEST  FOR  SIGN  CHANGE 

99  IF ( B ( K, 1 ) *B ( KM1 ,  1 ) . GE . 0 )  GO  TO  7 

94  NR-NR+ 1 
7  CONTINUE 

C  IF  ALL  ROOTS  ARE  TO  THE  RIGHT  AND  CONVERGENCE  HAS  NOT 
C  STARTED,  SHIFT  AXIS  LEFT  BY  A NO 
IF(N.LE.NR)  GO  TO  24 

21  IF  (KS.GT.O)  GO  TO  26 

22  AL=ALf ANO 
GO  TO  44 

24  IF(KS.NE.O)  GO  TO  26 

25  NO=NR 

26  IF(NO.LE.O)  GO  TO  1 
31  IF(NO.GT.NR)  GO  TO  34 
33  IF(KK.NE.O)  GO  TO  35 
37  AL=AL+ALO 

DA=AL0/2 . DO 


ee  ot  oo  (:a.ao.((i .x)a)8aAa)^i 
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GO  TO  38 
35  AL=AL+DA 
41  DA=DA/2.D0 
GO  TO  38 
34  AL=AL-DA 
KK=1 


IF(DABS(DA/AL) 

.LE.E3) 

GO 

C\J 

^r 

0 

E-i 

43 

IF (DABS (DA ) . GT 

■  El) 

GO 

TO 

41 

IF (DABS (DA ) . LE 

•  El) 

GO 

TO 

42 

38 

KS=  1 

CALCULATE  AA ( 2 ) . . 

.  .  . AA(Ntl) 

44 

DO  30  1=2, N1 

AA(l)=A(l) 

N2I=N2-I 
IM1=I-1 
DO  30  J=1,IM1 
IMJ=I- J 
N1J=N1-J 

30  A A ( I) =AA ( I) f C ( N1 J, N2I )  *A  ( J ) *( AL**IMJ ) 

GO  TO  40 
42  NRR=NO~NR 
KK-0 

WRITE  (6,84) NRR , A L 

84  FORMAT  (1H0,I3,23H  ROOTS  WITH  REAL  PART=  D25.16/) 
L=N1-NRR 
LP=  ( L+-2  )/2 

C  TEST  FOR  0,2,4  OR  MORE  IMAGINARY  PARTS 
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IF(NRR.LE.l)  GO  TO  25 

51  IF(NRR.LE.3)  GO  TO  55 

52  IF(NRR. LE. 5)  GO  TO  53 

56  L2-L/2 
WRITE  (6,57) 

57  FORMAT  (48H  THE  IMAG  PARTS  ARE  THE  SORTS  OF  THE  ROOTS 
10F  AN/) 

WRITE (6,  58) L2 

58  FORMAT  (42H  AUXILIARY  EQTN  A ( I ) S**N. . . A (N+l )  OF  ORDER, 
113/) 

WRITE  (6,59)  (B(L,M),M=1,LP 

59  FORMAT  (4lH  WHOSE  COEFF  IN  DESCENDING  ORDER  OF  S  ARE,/ 
11X,  6D18. 10/) 

IF( L. LE. 2*L2)  GO  TO  25 

61  WRITE  (6,62) 

62  FORMAT  (26H  A  SINGLE  ROOT  ALSO  EXISTS/) 

GO  TO  25 

55  R1=DSQRT(DABS(B(L,2)/B(L,1))) 

R2=-R1 

WRITE  (6,85)  R1,R2 

85  FORMAT  (33H  2  OF  THESE  HAVE  COMPLEX  PARTS  AT/2D25.16/) 
GO  TO  25 

53  R1=-B(L,2)/(2.D0*B(L,1)) 

R3=DSQRT(R1*R1-B(L,3)/B(L,  1) ) 

R4=DSQRT(DABS(R1+R3) ) 

R2=DSQRT(DABS(R1-R3) ) 

Rl=-R2 

R3=-R4 


BTooff  a:  t  w  2THP2  sht  sha  stfah  oayiI  $i,:t  -tB-0  a*-.-  o  i 

/ 

HSOHO  r-iO  (I+M)A. .  .H**a(l)  MTP3  YHAIJIXUA  H2*)  1  tHO'I  85 

. 
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WRITE  (6,86)  R2,R1,R4,R3 

86  FORMAT  (33H  4  OF  THESE  HAVE  COMPLEX  PARTS  AT/4D25.16/) 
GO  TO  25 
END 


' 


CHAPTER  VI 


VALIDATION  OF  RESULTS  OBTAINED  USING  ROUTH'S  CRITERION  FOR 
POLYNOMIAL  ROOT  SOLUTION 

This  chapter  compares  the  results  obtained  using  the 
computer  program  of  Chapter  V  with  those  obtained  using  a 
program  from  the  SHARE  Library  of  the  Computing  Center  at 
the  University  of  Alberta. 

6.1  SHARE  LIBRARY  PROGRAM  C2-3332 

This  program  can  be  used  to  find  the  roots  of  a  given 
polynomial  of  any  order  or  to  generate  a  polynomial  from  the 
given  roots.  When  solving  for  the  roots  of  a  polynomial,  the 
program  generates  a  polynomial  from  the  solved  roots  which 
can  be  compared  to  the  given  polynomial.  This  serves  as  a 
built-in  accuracy  check  of  the  calculated  roots. 

The  roots  are  found  using  a  combination  of  Newton  and 
Muller  approximations.  Muller  iterations  provide  a  first 
approximation  to  the  root  and  the  Newton  iterations  continue 
the  solution  until  the  root  is  found  to  the  desired  accuracy. 
The  root,  roots  if  complex  conjugates,  is  extracted  from  the 
polynomial  and  the  solution  continues  with  the  reduced 
polynomial. 

All  calculations  are  carried  out  using  double 
precision  arithmetic,  but  the  solutions  are  printed  out  in 
single  precision  arithmetic. 

6.2  RESULTS  OBTAINED  USING  ROUTH'S  CRITERION  PROGRAM 


The  SHARE  Library  program  C2-3332  was  used  as  a 
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standard  to  check  the  accuracy  of  the  program  based  on  Routh's 
criterion . 

An  ill-conditioned  polynomial  was  chosen  to  test  the 
accuracy  of  the  programs,  as  it  is  one  of  the  most  difficult 
types  of  polynomials  to  solve.  Small  changes  in  the  coeffic¬ 
ients  of  an  ill-conditioned  polynomial  will  cause  large 
changes  in  the  values  of  the  roots.  The  ill-conditioned 
polynomial  shown  below  was  chosen  from  the  article  "The  Eval¬ 
uation  of  the  Zeros  of  Ill-Conditioned  Polynomials"  by  J.  H. 
Wilkinson  printed  in  Numerische  Ma thematik,  Part  1,  Springer- 
Verlag,  1959- 

2.03253121s16  f  3.4356o48s16  *  25. 1783048s12*  +  37.651096s18  + 
128.218748s12  f  166.44768s11  4-  345.07256s10  f  378.908s9  f 
524.327s8  t  468.88s7  f  443.576s6  f  304. 08s6  4-190. 68s21  4- 
89.6s8  4-  32.8s2  4-  8s  4-  1  =  0  . (6-1) 

The  solutions  given  in  the  literature  were  correct 

to  15  decimal  places,  but  were  rounded-off  to  9  decimals  as 

shown  in  Table  6.1.  Solutions  of  the  polynomial  using  the 

SHARE  program  and  the  Routh’s  criterion  program  are  also  shown 

in  Table  6.1.  The  values  used  in  the  Routh's  criterion 

—  8  — 10 

program  for  El  and  E3  were  1x10  and  1x10  respectively. 

Table  6.1  shows  that  the  roots  obtained  using  the 
SHARE  program  are  in  agreement  with  the  correct  values  for 
the  roots  to  eight  significant  figures.  The  results  obtained 
using  the  Routh  program  were  in  agreement  to  nine  significant 
figures  with  the  correct  values,  except  for  roots  5  and  6 
which  only  agreed  to  eight  significant  figures,  but  only  14 
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ROOTS 

CORRECT  VALUE 

OF  ROOT 

RESULTS  USING 
SHARE  PROGRAM 
C2-3332 

RESULTS  USING 

ROUTH  CRITERION 
PROGRAM 

1,2 

-0.293504529 

-jo. 143499296 

-0.29350453 

±J0. 14349930 

-0.293504529 

-JO. 143499296 

3,4 

-0.224470057 

ijo. 450927958 

-0.22447006 

^0.45092796 

-0.224470057 

t jo. 450927958 

5,6 

-0.147623780 

-JO. 771757201 

-0.14762378 

-  JO. 77175720 

-0.147623780 

4  JO. 771757200 

7,8 

-0.0900399887 

tji. 06119206 

-0.090039989 

ijl. 0611921 

-0.0900399887 

ijl. 06119206 

9,10 

-0.0508644356 

±jl. 29691128 

-0.050864436 

*J1. 2969113 

-0.0508644356 

- jl. 29691128 

11  j  12 

-0.0256687105 

-3 1.47437714 

-0.025668711 

-Jl. 4743771 

-0.0256687105 

4 Jl. 47437714 

13,14 

-0.0104935501 

tj 1.59629550 

-0.010493550 

-Jl. 5962955 

-0.0104935501 

4  J  1.59629550 

15,16 

-0.00248920244 

tj 1.66712036 

-0.0024892024 

1.6671204 

PROGRAM 

TERMINATED 

TABLE  6.1  SOLUTIONS  FOR  ILL-CONDITIONED  POLYNOMIAL  (6-1) 

of  the  16  roots  were  calculated.  An  excessive  number  of 
floating  point  traps  were  encountered  when  searching  for  roots 
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15  and  16  and  tine  calculations  were  thus  terminated.  Float¬ 
ing  point  traps  resulting  in  the  output  of  underflow  messages 
were  also  encountered  when  the  computer  was  searching  for 
roots  11,  12,  13  and  14,  but  there  were  not  enough  of  them 
to  terminate  the  calculations. 

To  further  compare  the  results  of  the  two  programs, 
a  second  polynomial  shown  below,  not  ill-conditioned,  was 
solved  and  the  roots  obtained  listed  in  Table  6.2. 

s9  +■  21.077365s8  +  173.21313s7  *  758.54868s8  +■ 
2185.2366s5  f  4804.673s21  +  7758.6178s3  + 

8765.4409s2  +  7417.3856s  +  1692.535  =  0  . (6-2) 


ROOTS 

RESULTS  USING 
SHARE  PROGRAM 
C2-3332 

RESULTS  USING 
ROTJTH  CRITERION 
PROGRAM 

1 

-7.7770889 

-7.7770889 

2 

-5. 4321411 

-5.4321411 

3,4 

-3.2139902 

iJO. 12334476 

-3.2139902 

■t  JO.  12334476 

5,6 

-0.49999812 

-jl. 7321397 

-0.49999812 

iji. 7321397 

7 

-0.32167519 

-0.32167519 

8,9 

-0.059241578 

-Jl. 9236846 

-0.059241578 

ijl. 9236846 

TABLE  6.2  SOLUTIONS  FOR  POLYNOMIAL  (6-2) 


■ 
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The  values  of  the  roots  obtained  using  Routin' s 
criterion  were  rounded-off  to  eight  significant  figures. 

An  examination  of  Table  6.2  shows  that  the  results 
obtained  in  both  cases  are  identical.  Although  all  the  roots 
were  calculated  using  the  Routh's  criterion  program,  floating 
point  traps  were  encountered  while  the  computer  was  searching 
for  root  number  7. 

The  reason  for  the  floating  point  traps  is  the 
extremely  small  numbers  calculated  when  determining  the 
coefficients  AA(l)  of  the  polynomial  with  shifted  axis. 

These  occur  when  the  shifted  axis  is  close  to  the  original 
axis,  that  is  AL  very  small. 

For  the  case  of  the  ill-conditioned  polynomial,  when 

the  computer  was  searching  for  the  roots  with  real  part  at 

-0.00248920244,  the  axis  would  be  between  -0.001  and  -0.003 

sometime  during  the  search.  Raising  any  number  in  this  range 

-38 

to  the  sixteenth  power,  produces  a  number  less  than  1x10  ^ 
which  will  interrupt  the  computer.  The  floating  point  trap 
substitutes  a  zero  and  the  calculations  proceed,  but  the 
amount  of  axis  shift,  AL,  eventually  becomes  small  enough 
to  cause  another  underflow  and  thus  the  solution  is  terminated 
before  completion. 

In  the  case  of  the  polynomial  (6-2),  the  underflows 
were  also  caused  by  numbers  less  than  1x10“^®,  They  arose 
when  the  axis  was  shifted  0.5  units  to  the  right,  AL0r0.5^ 
from  -0.49999812  thus  setting  AL  at  0. 00000188.  Raising  this 
number  to  the  ninth  power  caused  the  underflow.  But  in  this 
incidence  the  floating  point  trap  substituted  zero  and  the 
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calculations  proceeded  until  the  solution  was  completed. 

6.3  TIME  REQUIRED  FOR  CALCULATIONS 

The  time  required  for  calculations,  that  is  the  actual 
program  running  time  or  object  time,  and  the  total  time 
required  to  process  the  program  are  shown  in  Table  6.3  for 
both  the  Routh  program  and  the  SHARE  program. 


EQUATION 

ROUTH  PROGRAM 

SHARE  PROGRAM 

OBJECT 

TIME 

TOTAL 

TIMS 

OBJECT 

TIME 

TOTAL 

TIME 

6-1 

27 

83 

11 

84 

6-2 

22 

58 

4 

52 

TABLE  6.3  TIME  (IN  SECONDS)  REQUIRED  FOR  CALCULATIONS 

Although  the  object  time  of  the  SHARE  program  is  much 
less  than  that  of  the  Routh  program  for  the  polynomial 
solutions  considered,  the  total  times  required  to  run  the 
programs  are  approximately  the  same  in  both  cases.  Thus,  if 
time  was  a  factor  in  choosing  a  program  for  a  particular 
solution,  the  SHARE  program  offers  the  advantage  of  a  shorter 
computer  running  time. 

6.4  CONCLUSIONS 


The  program  based  on  Routh* s  stability  criterion  will 
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solve  for  the  roots  of  almost  any  polynomial  to  the  accuracy 
stipulated  in  the  input  data.  Polynomials  for  which  all  the 
roots  may  not  be  calculated  are  those  which  do  not  satisfy 
the  conditions  (AL)n  >  10”^  and  (AL)n  <  10^  where  AL  is  the 
real  part  of  the  root  and  N  is  the  degree  of  the  polynomial. 
Failure  of  the  roots  to  satisfy  these  conditions  causes  an 
excessive  number  of  underflows  or  overflows  respectively 
and  the  calculations  are  terminated  in  either  case. 

In  engineering  practice,  polynomials  with  roots  not 
satisfying  the  above  conditions  are  not  normally  encountered. 
Thus  the  engineer  can  use  this  program  with  confidence  that 
the  solutions  will  be  correct  and  complete.  If  by  chance 
any  roots  did  lie  outside  these  extremes,  the  roots  calculated 
and  printed  out  would  be  correct,  but  the  solution  would  not 
be  complete. 

This  program  may  appeal  to  the  electrical  engineer 
because  it  is  based  on  Routin' s  stability  criterion  with 
which  he  is  familiar.  This  allows  him  to  vary  the  input  par- 
amaters  with  an  understanding  of  the  effects  the  changes 
will  have  on  the  calculations. 
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APPENDIX  A 


FLOWCHART  FOR  PROGRAM  OF  CHAPTER  III 
A . 1  MAIN  PROGRAM 


(  eqns(a,b,e,f,n,t7~) 


f  N,  T 

N 

CALL 


ARRAYS 

N,T 


T 


F(5,I)=F(1,6) 


3 

J=2,N 
I 


F(5, J)=F(J,5) 


4 

1=1,  N 


F(I*5,1)SF.(I,1)*F(5,1) 


F(6,2)=F(1,6) 
5 

1=2,  N 

i 


p(i+5,2)=p(i,5) 


> - > - 

A 

F(6,3)=F(1,2)*F(5,1)<-f(1,3)*F(5,2)  + 
f(i,A)*f(5,3)+p(i,5)*f(5,^) 
Zd 

1=2,  N 


B 


Ill  AHO  TO  MA8G0  1  80'"  IA  *')  '  I' 
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^ _ 

F(I+5,3)=F 
+  F 

fcli 

*F(5, 2) 
*f(5,4) 

+F(I,3)*F(5,3) 

F(I+9,  J)=F(I,l)*F(6,  J) 


11 
1=1,3 


F(lO,I+3)=F(l,2)*F(6,l)+F(l,3)*F(7,I 

+F(1,4)*F(8,I)+f(1,5)*F(9,I 


F(10,7)=F(1,6) 
T 


F(l+10, J+3)=F(I+1,2)*F(7, J)+F(l+1,3)*F(8,  J) 
+F( 1+1,4 ) *F(9, J) 


F(I  +  10,7)=F(I+1,5)| 


F(I+13, J)=F(l,l)*F(lO, J) 


19 

1=1,4 


F(I+13,8)-F(I+9,7) 


■ 
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48 

1=1,3 

J=l,  8 
.  1 


A(I,  j)=B(l+10,  J+8) 


Kp) 

C  GO  TO  80  > 


E 


49K 


E(4,l 
E  5,1 
E(6,l 


=E(1,5 
=E  2,4 
=E(3, 4 


50 

1=1,3 


E(l+6,l)=E(l,l)*E(4,l) 

So: 


52 

1=1,3 


E(  14-6 , 2 )  =e(i+3j  1 ) 
52) 


E(7,3)=E(l,2)»E(4,l)+E(l,3)*E(5,l)*E(l,4)»E(6,l) 


54  ' 

1=1,2 


E(l+7,3)=E(l+l,2)*E(5,l)+E(l+l,3)*E(6,l) 


56 

1=1,3 

J=l,3 


E(lf9,  J)=E(I,1)*E(7,  J) 


F 


' 
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58 

1=1,3 


E(10,I+3)=E(1,2)*E(7,I)+E(1,3)*E(8,I) 

_ +e(1,4)»e(9,I) _ 


60 

1=1,3 


E(I+9,7)=E(3>3,1) 


62 

1=1,2 

J=l,3 


E(l+10, J+3)=E(I+1,2)*E(8, J)+E(l+1,3)*E(9, J) 

>62, 


66 

1=1, N 

J=l,7 


I 


A(I, J)=E(l+9,  J) 


E 


82 

1=1,2 

J=l,  8 


A(l+3,  J)=A(I+1,  J)*A(1,7) 


84 

1=2,3 

J=l,8 


I 


A ( 1+4 ,  j)=A(l,j)*A(l,7) 

^4) 


G 
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A(l+7, J)=A(I+3, J)-A(l+5, J) 


87 

J-1,8 


A(10,J)  =  A(8,  J)*A(9>5) 


88 

J-1,8 


A(ll, J)=A(9, J)*A(8,5) 


89 

J=l,8 

i 


A(12, J)=A(10, J)-A(ll, J) 


CK0=A( 12, 8) /A  (12, 6) 

CK1= ( A (8, 8) -A (8,  6) *CK0)/(  A(8, 4 ) *CK0- A (8, 5 ) ) 
CNK2=A(1,8)-A(1,4)*CK0*CK1-A(1,5)*CK1-A(1,6)*CK0 
CDK2=A(1,1)*CK0*CK1+A(1,2)*CK1+A(1,3)*CK0+A(1,7) 
CK2=CNK2/CDK2 _ 


T 


CK01=CK0*CK1 

CK02=CK0*CK2 

CK12=CK1*CK2 

CK012=CK0*CK1*CK2 

CKN3=B( 1, 16) -B( 1,9 ) *CK012-B( 1, 10) *CK12-B( 1, 11 ) *CK02- 
B(l,12)*CK01-B(l,13)*CKl-B(l,l4)*CK0-B(l,15)*CK2 
CKD3=B( 1,1) *CK012+B( 1 . 2 ) *CK12+B (1,3) *CK02+B( 1,4) *CK01+ 
B(1,5)*CK1+B(1,6)*CK04-B(1,7)*CK2+B(1,8) 
CK3-CKN3/CKD3 


■ 
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HO=CKO*ERO 

H1=CK1*ER1 

H2=CK2*ER2 


ER0,ER1,ER2 
HO,  HI ,  H2 _ 


(  call  exit"") 


ERO=B! 

ER1=B( 

[1, 16] 
1,16 

-f(5> i) *cko 

ER2=BI 

1, 16 

-f(6, 1) *cko*cki-f(6, 2) *cki-f(6,  3)*CK0 

ER3=BI 

1,16 

-F( 10, 1 ) *CK012-Ff 10, 2) *CK12-F( 10,  3 ) *CK02- 

F( 10, 4 ) *CK01-F( 10, 5 ) *CK1-F( 10, 6) *CKO-F( 10, 7) *CK2 

HO=CKO*ERO 

H1=CK1*ER1 

H2=CK2*ER2 

H3ZCK3*ER3 


ERO,ERl,ER2,ER3 
HO, HI, H2, H3 _ 


a=SHH| 

a=era 


* 


-  72 


A. 2  SUBROUTINE  ARRAYS 


(  SUBROUTINE  ARRAYS  (N,T) 


PUNCH  IN  THE  APPROPRIATE 
COEFFICIENTS  AND  FINAL 
CONDITIONS  AS  OUTLINED 
IN  THE  PROGRAM  IN 
CHAPTER  III 


(  RETURN  ) 


END 


APPENDIX  B 


FLOW  CHART  FOR  PROGRAM  TO  FIND  THE  ROOTS  OF  A  POLYNOMIAL 
USING  ROUTH'S  STABILITY  CRITERION 
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76  - 


-  77  - 


78 
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T 


ALr 

ALfDA 

© 

DA=DA/2. 

—  - 1 - 

(  GO 

TO 

38 

@ 

rv  ■  - 

AL= 

KK= 

AL 

1 

-DA 

DA 

TF  T7TV— 

T 

AL 

•  ±-iHi .  ro 

AA  ( I ) =A ( I ) 
N2I=N2-I 
IM1=I-1 

30 

J= 1 , IMl 


M 


T 


IMJ=I-J 
N1  J=N1-  J 

AA(I)=AA(I)+C(N1J,N2I)*A( J)*(AL**IMJ) 


. 
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A  SINGLE  ROOT  ALSO 
EXISTS 


(“  GO  TO  25  ) 


Rl=DSQRT(DABS(B(L,2)/B(L,l) ) ) 
R2=-R1 


R1,R2 

2  OF  THESE  HAVE  COMPLEX 
PARTS  AT  _ 


(  GO  TC 

in 

OJ 

■  / - N. 

£ 

"3^  < 

r 

R1=-B(L,2)/(2*B(L,1)) 
R3=DSQRT(R1*R1-B(L,3  )/b(l, 1) ) 
R^=DSQRT  DABS(R1+R3) 
R2=DSQRT(DABS(R1-R3) ) 

Rl=-R2 

R3=-R'4 

R2,R1,F 

!4,R3 

T 


OF  THESE  HAVE  COMPLEX 
PARTS  AT 


(  GO  TO  28  ) 
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